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
  45 /*
  46  * This file includes the code, contributed by Simon Perreault
  47  * (the function icvMedianBlur_8u_O1)
  48  *
  49  * Constant-time median filtering -- http://nomis80.org/ctmf.html
  50  * Copyright (C) 2006 Simon Perreault
  51  *
  52  * Contact:
  53  *  Laboratoire de vision et systemes numeriques
  54  *  Pavillon Adrien-Pouliot
  55  *  Universite Laval
  56  *  Sainte-Foy, Quebec, Canada
  57  *  G1K 7P4
  58  *
  59  *  perreaul@gel.ulaval.ca
  60  */
  61
  62 namespace cv
  63 {
  64
  65 /****************************************************************************************\
  66                                          Box Filter
  67 \****************************************************************************************/
  68
  69 template<typename T, typename ST> struct RowSum : public BaseRowFilter
  70 {
  71     RowSum( int _ksize, int _anchor )
  72     {
  73         ksize = _ksize;
  74         anchor = _anchor;
  75     }
  76
  77     void operator()(const uchar* src, uchar* dst, int width, int cn)
  78     {
  79         const T* S = (const T*)src;
  80         ST* D = (ST*)dst;
  81         int i = 0, k, ksz_cn = ksize*cn;
  82
  83         width = (width - 1)*cn;
  84         for( k = 0; k < cn; k++, S++, D++ )
  85         {
  86             ST s = 0;
  87             for( i = 0; i < ksz_cn; i += cn )
  88                 s += S[i];
  89             D[0] = s;
  90             for( i = 0; i < width; i += cn )
  91             {
  92                 s += S[i + ksz_cn] - S[i];
  93                 D[i+cn] = s;
  94             }
  95         }
  96     }
  97 };
  98
  99
 100 template<typename ST, typename T> struct ColumnSum : public BaseColumnFilter
 101 {
 102     ColumnSum( int _ksize, int _anchor, double _scale )
 103     {
 104         ksize = _ksize;
 105         anchor = _anchor;
 106         scale = _scale;
 107         sumCount = 0;
 108     }
 109
 110     void reset() { sumCount = 0; }
 111
 112     void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 113     {
 114         int i;
 115         ST* SUM;
 116         bool haveScale = scale != 1;
 117         double _scale = scale;
 118
 119         if( width != (int)sum.size() )
 120         {
 121             sum.resize(width);
 122             sumCount = 0;
 123         }
 124
 125         SUM = &sum[0];
 126         if( sumCount == 0 )
 127         {
 128             for( i = 0; i < width; i++ )
 129                 SUM[i] = 0;
 130             for( ; sumCount < ksize - 1; sumCount++, src++ )
 131             {
 132                 const ST* Sp = (const ST*)src[0];
 133                 for( i = 0; i <= width - 2; i += 2 )
 134                 {
 135                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 136                     SUM[i] = s0; SUM[i+1] = s1;
 137                 }
 138
 139                 for( ; i < width; i++ )
 140                     SUM[i] += Sp[i];
 141             }
 142         }
 143         else
 144         {
 145             CV_Assert( sumCount == ksize-1 );
 146             src += ksize-1;
 147         }
 148
 149         for( ; count--; src++ )
 150         {
 151             const ST* Sp = (const ST*)src[0];
 152             const ST* Sm = (const ST*)src[1-ksize];
 153             T* D = (T*)dst;
 154             if( haveScale )
 155             {
 156                 for( i = 0; i <= width - 2; i += 2 )
 157                 {
 158                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 159                     D[i] = saturate_cast<T>(s0*_scale);
 160                     D[i+1] = saturate_cast<T>(s1*_scale);
 161                     s0 -= Sm[i]; s1 -= Sm[i+1];
 162                     SUM[i] = s0; SUM[i+1] = s1;
 163                 }
 164
 165                 for( ; i < width; i++ )
 166                 {
 167                     ST s0 = SUM[i] + Sp[i];
 168                     D[i] = saturate_cast<T>(s0*_scale);
 169                     SUM[i] = s0 - Sm[i];
 170                 }
 171             }
 172             else
 173             {
 174                 for( i = 0; i <= width - 2; i += 2 )
 175                 {
 176                     ST s0 = SUM[i] + Sp[i], s1 = SUM[i+1] + Sp[i+1];
 177                     D[i] = saturate_cast<T>(s0);
 178                     D[i+1] = saturate_cast<T>(s1);
 179                     s0 -= Sm[i]; s1 -= Sm[i+1];
 180                     SUM[i] = s0; SUM[i+1] = s1;
 181                 }
 182
 183                 for( ; i < width; i++ )
 184                 {
 185                     ST s0 = SUM[i] + Sp[i];
 186                     D[i] = saturate_cast<T>(s0);
 187                     SUM[i] = s0 - Sm[i];
 188                 }
 189             }
 190             dst += dststep;
 191         }
 192     }
 193
 194     double scale;
 195     int sumCount;
 196     vector<ST> sum;
 197 };
 198
 199
 200 template<> struct ColumnSum<int, uchar> : public BaseColumnFilter
 201 {
 202     ColumnSum( int _ksize, int _anchor, double _scale )
 203     {
 204         ksize = _ksize;
 205         anchor = _anchor;
 206         scale = _scale;
 207         sumCount = 0;
 208     }
 209
 210     void reset() { sumCount = 0; }
 211
 212     void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 213     {
 214         int i;
 215         int* SUM;
 216         bool haveScale = scale != 1;
 217         double _scale = scale;
 218
 219         #if CV_SSE2
 220             bool haveSSE2 = checkHardwareSupport(CV_CPU_SSE2);
 221         #endif
 222
 223         if( width != (int)sum.size() )
 224         {
 225             sum.resize(width);
 226             sumCount = 0;
 227         }
 228
 229         SUM = &sum[0];
 230         if( sumCount == 0 )
 231         {
 232             memset((void*)SUM, 0, width*sizeof(int));
 233             for( ; sumCount < ksize - 1; sumCount++, src++ )
 234             {
 235                 const int* Sp = (const int*)src[0];
 236                 i = 0;
 237                 #if CV_SSE2
 238                 if(haveSSE2)
 239                 {
 240                     for( ; i < width-4; i+=4 )
 241                     {
 242                         __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
 243                         __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
 244                         _mm_storeu_si128((__m128i*)(SUM+i),_mm_add_epi32(_sum, _sp));
 245                     }
 246                 }
 247                 #endif
 248                 for( ; i < width; i++ )
 249                     SUM[i] += Sp[i];
 250             }
 251         }
 252         else
 253         {
 254             CV_Assert( sumCount == ksize-1 );
 255             src += ksize-1;
 256         }
 257
 258         for( ; count--; src++ )
 259         {
 260             const int* Sp = (const int*)src[0];
 261             const int* Sm = (const int*)src[1-ksize];
 262             uchar* D = (uchar*)dst;
 263             if( haveScale )
 264             {
 265                 i = 0;
 266                 #if CV_SSE2
 267                 if(haveSSE2)
 268                 {
 269                     const __m128 scale4 = _mm_set1_ps((float)_scale);
 270                     for( ; i < width-8; i+=8 )
 271                     {
 272                         __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
 273                         __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));
 274
 275                         __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 276                                                      _mm_loadu_si128((const __m128i*)(Sp+i)));
 277                         __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
 278                                                       _mm_loadu_si128((const __m128i*)(Sp+i+4)));
 279
 280                         __m128i _s0T = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
 281                         __m128i _s0T1 = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s01)));
 282
 283                         _s0T = _mm_packs_epi32(_s0T, _s0T1);
 284
 285                         _mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));
 286
 287                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 288                         _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
 289                     }
 290                 }
 291                 #endif
 292                 for( ; i < width; i++ )
 293                 {
 294                     int s0 = SUM[i] + Sp[i];
 295                     D[i] = saturate_cast<uchar>(s0*_scale);
 296                     SUM[i] = s0 - Sm[i];
 297                 }
 298             }
 299             else
 300             {
 301                 i = 0;
 302                 #if CV_SSE2
 303                 if(haveSSE2)
 304                 {
 305                     for( ; i < width-8; i+=8 )
 306                     {
 307                         __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
 308                         __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));
 309
 310                         __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 311                                                      _mm_loadu_si128((const __m128i*)(Sp+i)));
 312                         __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
 313                                                       _mm_loadu_si128((const __m128i*)(Sp+i+4)));
 314
 315                         __m128i _s0T = _mm_packs_epi32(_s0, _s01);
 316
 317                         _mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));
 318
 319                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 320                         _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
 321                     }
 322                 }
 323                 #endif
 324
 325                 for( ; i < width; i++ )
 326                 {
 327                     int s0 = SUM[i] + Sp[i];
 328                     D[i] = saturate_cast<uchar>(s0);
 329                     SUM[i] = s0 - Sm[i];
 330                 }
 331             }
 332             dst += dststep;
 333         }
 334     }
 335
 336     double scale;
 337     int sumCount;
 338     vector<int> sum;
 339 };
 340
 341 template<> struct ColumnSum<int, short> : public BaseColumnFilter
 342 {
 343     ColumnSum( int _ksize, int _anchor, double _scale )
 344     {
 345         ksize = _ksize;
 346         anchor = _anchor;
 347         scale = _scale;
 348         sumCount = 0;
 349     }
 350
 351     void reset() { sumCount = 0; }
 352
 353     void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 354     {
 355         int i;
 356         int* SUM;
 357         bool haveScale = scale != 1;
 358         double _scale = scale;
 359
 360         #if CV_SSE2
 361             bool haveSSE2 = checkHardwareSupport(CV_CPU_SSE2);
 362         #endif
 363
 364         if( width != (int)sum.size() )
 365         {
 366             sum.resize(width);
 367             sumCount = 0;
 368         }
 369         SUM = &sum[0];
 370         if( sumCount == 0 )
 371         {
 372             memset((void*)SUM, 0, width*sizeof(int));
 373             for( ; sumCount < ksize - 1; sumCount++, src++ )
 374             {
 375                 const int* Sp = (const int*)src[0];
 376                 i = 0;
 377                 #if CV_SSE2
 378                 if(haveSSE2)
 379                 {
 380                     for( ; i < width-4; i+=4 )
 381                     {
 382                         __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
 383                         __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
 384                         _mm_storeu_si128((__m128i*)(SUM+i),_mm_add_epi32(_sum, _sp));
 385                     }
 386                 }
 387                 #endif
 388                 for( ; i < width; i++ )
 389                     SUM[i] += Sp[i];
 390             }
 391         }
 392         else
 393         {
 394             CV_Assert( sumCount == ksize-1 );
 395             src += ksize-1;
 396         }
 397
 398         for( ; count--; src++ )
 399         {
 400             const int* Sp = (const int*)src[0];
 401             const int* Sm = (const int*)src[1-ksize];
 402             short* D = (short*)dst;
 403             if( haveScale )
 404             {
 405                 i = 0;
 406                 #if CV_SSE2
 407                 if(haveSSE2)
 408                 {
 409                     const __m128 scale4 = _mm_set1_ps((float)_scale);
 410                     for( ; i < width-8; i+=8 )
 411                     {
 412                         __m128i _sm   = _mm_loadu_si128((const __m128i*)(Sm+i));
 413                         __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));
 414
 415                         __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 416                                                      _mm_loadu_si128((const __m128i*)(Sp+i)));
 417                         __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
 418                                                       _mm_loadu_si128((const __m128i*)(Sp+i+4)));
 419
 420                         __m128i _s0T  = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
 421                         __m128i _s0T1 = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s01)));
 422
 423                         _mm_storeu_si128((__m128i*)(D+i), _mm_packs_epi32(_s0T, _s0T1));
 424
 425                         _mm_storeu_si128((__m128i*)(SUM+i),_mm_sub_epi32(_s0,_sm));
 426                         _mm_storeu_si128((__m128i*)(SUM+i+4), _mm_sub_epi32(_s01,_sm1));
 427                     }
 428                 }
 429                 #endif
 430                 for( ; i < width; i++ )
 431                 {
 432                     int s0 = SUM[i] + Sp[i];
 433                     D[i] = saturate_cast<short>(s0*_scale);
 434                     SUM[i] = s0 - Sm[i];
 435                 }
 436             }
 437             else
 438             {
 439                 i = 0;
 440                 #if CV_SSE2
 441                 if(haveSSE2)
 442                 {
 443                     for( ; i < width-8; i+=8 )
 444                     {
 445
 446                         __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
 447                         __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));
 448
 449                         __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 450                                                      _mm_loadu_si128((const __m128i*)(Sp+i)));
 451                         __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
 452                                                       _mm_loadu_si128((const __m128i*)(Sp+i+4)));
 453
 454                         _mm_storeu_si128((__m128i*)(D+i), _mm_packs_epi32(_s0, _s01));
 455
 456                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 457                         _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
 458                     }
 459                 }
 460                 #endif
 461
 462                 for( ; i < width; i++ )
 463                 {
 464                     int s0 = SUM[i] + Sp[i];
 465                     D[i] = saturate_cast<short>(s0);
 466                     SUM[i] = s0 - Sm[i];
 467                 }
 468             }
 469             dst += dststep;
 470         }
 471     }
 472
 473     double scale;
 474     int sumCount;
 475     vector<int> sum;
 476 };
 477
 478
 479 template<> struct ColumnSum<int, ushort> : public BaseColumnFilter
 480 {
 481     ColumnSum( int _ksize, int _anchor, double _scale )
 482     {
 483         ksize = _ksize;
 484         anchor = _anchor;
 485         scale = _scale;
 486         sumCount = 0;
 487     }
 488
 489     void reset() { sumCount = 0; }
 490
 491     void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
 492     {
 493         int i;
 494         int* SUM;
 495         bool haveScale = scale != 1;
 496         double _scale = scale;
 497         #if CV_SSE2
 498                 bool haveSSE2 =  checkHardwareSupport(CV_CPU_SSE2);
 499         #endif
 500
 501         if( width != (int)sum.size() )
 502         {
 503             sum.resize(width);
 504             sumCount = 0;
 505         }
 506         SUM = &sum[0];
 507         if( sumCount == 0 )
 508         {
 509             memset((void*)SUM, 0, width*sizeof(int));
 510             for( ; sumCount < ksize - 1; sumCount++, src++ )
 511             {
 512                 const int* Sp = (const int*)src[0];
 513                 i = 0;
 514                 #if CV_SSE2
 515                 if(haveSSE2)
 516                 {
 517                     for( ; i < width-4; i+=4 )
 518                     {
 519                         __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
 520                         __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
 521                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_add_epi32(_sum, _sp));
 522                     }
 523                 }
 524                 #endif
 525                 for( ; i < width; i++ )
 526                     SUM[i] += Sp[i];
 527             }
 528         }
 529         else
 530         {
 531             CV_Assert( sumCount == ksize-1 );
 532             src += ksize-1;
 533         }
 534
 535         for( ; count--; src++ )
 536         {
 537             const int* Sp = (const int*)src[0];
 538             const int* Sm = (const int*)src[1-ksize];
 539             ushort* D = (ushort*)dst;
 540             if( haveScale )
 541             {
 542                 i = 0;
 543                 #if CV_SSE2
 544                 if(haveSSE2)
 545                 {
 546                     const __m128 scale4 = _mm_set1_ps((float)_scale);
 547                     const __m128i delta0 = _mm_set1_epi32(0x8000);
 548                     const __m128i delta1 = _mm_set1_epi32(0x80008000);
 549
 550                     for( ; i < width-4; i+=4)
 551                     {
 552                         __m128i _sm   = _mm_loadu_si128((const __m128i*)(Sm+i));
 553                         __m128i _s0   = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 554                                                       _mm_loadu_si128((const __m128i*)(Sp+i)));
 555
 556                         __m128i _res = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
 557
 558                         _res = _mm_sub_epi32(_res, delta0);
 559                         _res = _mm_add_epi16(_mm_packs_epi32(_res, _res), delta1);
 560
 561                         _mm_storel_epi64((__m128i*)(D+i), _res);
 562                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 563                     }
 564                 }
 565                 #endif
 566                 for( ; i < width; i++ )
 567                 {
 568                     int s0 = SUM[i] + Sp[i];
 569                     D[i] = saturate_cast<ushort>(s0*_scale);
 570                     SUM[i] = s0 - Sm[i];
 571                 }
 572             }
 573             else
 574             {
 575                 i = 0;
 576                 #if  CV_SSE2
 577                 if(haveSSE2)
 578                 {
 579                     const __m128i delta0 = _mm_set1_epi32(0x8000);
 580                     const __m128i delta1 = _mm_set1_epi32(0x80008000);
 581
 582                     for( ; i < width-4; i+=4 )
 583                     {
 584                         __m128i _sm   = _mm_loadu_si128((const __m128i*)(Sm+i));
 585                         __m128i _s0   = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
 586                                                       _mm_loadu_si128((const __m128i*)(Sp+i)));
 587
 588                         __m128i _res = _mm_sub_epi32(_s0, delta0);
 589                         _res = _mm_add_epi16(_mm_packs_epi32(_res, _res), delta1);
 590
 591                         _mm_storel_epi64((__m128i*)(D+i), _res);
 592                         _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
 593                     }
 594                 }
 595                 #endif
 596
 597                 for( ; i < width; i++ )
 598                 {
 599                     int s0 = SUM[i] + Sp[i];
 600                     D[i] = saturate_cast<ushort>(s0);
 601                     SUM[i] = s0 - Sm[i];
 602                 }
 603             }
 604             dst += dststep;
 605         }
 606     }
 607
 608     double scale;
 609     int sumCount;
 610     vector<int> sum;
 611 };
 612
 613
 614 }
 615
 616 cv::Ptr<cv::BaseRowFilter> cv::getRowSumFilter(int srcType, int sumType, int ksize, int anchor)
 617 {
 618     int sdepth = CV_MAT_DEPTH(srcType), ddepth = CV_MAT_DEPTH(sumType);
 619     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(srcType) );
 620
 621     if( anchor < 0 )
 622         anchor = ksize/2;
 623
 624     if( sdepth == CV_8U && ddepth == CV_32S )
 625         return Ptr<BaseRowFilter>(new RowSum<uchar, int>(ksize, anchor));
 626     if( sdepth == CV_8U && ddepth == CV_64F )
 627         return Ptr<BaseRowFilter>(new RowSum<uchar, double>(ksize, anchor));
 628     if( sdepth == CV_16U && ddepth == CV_32S )
 629         return Ptr<BaseRowFilter>(new RowSum<ushort, int>(ksize, anchor));
 630     if( sdepth == CV_16U && ddepth == CV_64F )
 631         return Ptr<BaseRowFilter>(new RowSum<ushort, double>(ksize, anchor));
 632     if( sdepth == CV_16S && ddepth == CV_32S )
 633         return Ptr<BaseRowFilter>(new RowSum<short, int>(ksize, anchor));
 634     if( sdepth == CV_32S && ddepth == CV_32S )
 635         return Ptr<BaseRowFilter>(new RowSum<int, int>(ksize, anchor));
 636     if( sdepth == CV_16S && ddepth == CV_64F )
 637         return Ptr<BaseRowFilter>(new RowSum<short, double>(ksize, anchor));
 638     if( sdepth == CV_32F && ddepth == CV_64F )
 639         return Ptr<BaseRowFilter>(new RowSum<float, double>(ksize, anchor));
 640     if( sdepth == CV_64F && ddepth == CV_64F )
 641         return Ptr<BaseRowFilter>(new RowSum<double, double>(ksize, anchor));
 642
 643     CV_Error_( CV_StsNotImplemented,
 644         ("Unsupported combination of source format (=%d), and buffer format (=%d)",
 645         srcType, sumType));
 646
 647     return Ptr<BaseRowFilter>(0);
 648 }
 649
 650
 651 cv::Ptr<cv::BaseColumnFilter> cv::getColumnSumFilter(int sumType, int dstType, int ksize,
 652                                                      int anchor, double scale)
 653 {
 654     int sdepth = CV_MAT_DEPTH(sumType), ddepth = CV_MAT_DEPTH(dstType);
 655     CV_Assert( CV_MAT_CN(sumType) == CV_MAT_CN(dstType) );
 656
 657     if( anchor < 0 )
 658         anchor = ksize/2;
 659
 660     if( ddepth == CV_8U && sdepth == CV_32S )
 661         return Ptr<BaseColumnFilter>(new ColumnSum<int, uchar>(ksize, anchor, scale));
 662     if( ddepth == CV_8U && sdepth == CV_64F )
 663         return Ptr<BaseColumnFilter>(new ColumnSum<double, uchar>(ksize, anchor, scale));
 664     if( ddepth == CV_16U && sdepth == CV_32S )
 665         return Ptr<BaseColumnFilter>(new ColumnSum<int, ushort>(ksize, anchor, scale));
 666     if( ddepth == CV_16U && sdepth == CV_64F )
 667         return Ptr<BaseColumnFilter>(new ColumnSum<double, ushort>(ksize, anchor, scale));
 668     if( ddepth == CV_16S && sdepth == CV_32S )
 669         return Ptr<BaseColumnFilter>(new ColumnSum<int, short>(ksize, anchor, scale));
 670     if( ddepth == CV_16S && sdepth == CV_64F )
 671         return Ptr<BaseColumnFilter>(new ColumnSum<double, short>(ksize, anchor, scale));
 672     if( ddepth == CV_32S && sdepth == CV_32S )
 673         return Ptr<BaseColumnFilter>(new ColumnSum<int, int>(ksize, anchor, scale));
 674     if( ddepth == CV_32F && sdepth == CV_32S )
 675         return Ptr<BaseColumnFilter>(new ColumnSum<int, float>(ksize, anchor, scale));
 676     if( ddepth == CV_32F && sdepth == CV_64F )
 677         return Ptr<BaseColumnFilter>(new ColumnSum<double, float>(ksize, anchor, scale));
 678     if( ddepth == CV_64F && sdepth == CV_32S )
 679         return Ptr<BaseColumnFilter>(new ColumnSum<int, double>(ksize, anchor, scale));
 680     if( ddepth == CV_64F && sdepth == CV_64F )
 681         return Ptr<BaseColumnFilter>(new ColumnSum<double, double>(ksize, anchor, scale));
 682
 683     CV_Error_( CV_StsNotImplemented,
 684         ("Unsupported combination of sum format (=%d), and destination format (=%d)",
 685         sumType, dstType));
 686
 687     return Ptr<BaseColumnFilter>(0);
 688 }
 689
 690
 691 cv::Ptr<cv::FilterEngine> cv::createBoxFilter( int srcType, int dstType, Size ksize,
 692                     Point anchor, bool normalize, int borderType )
 693 {
 694     int sdepth = CV_MAT_DEPTH(srcType);
 695     int cn = CV_MAT_CN(srcType), sumType = CV_64F;
 696     if( sdepth <= CV_32S && (!normalize ||
 697         ksize.width*ksize.height <= (sdepth == CV_8U ? (1<<23) :
 698             sdepth == CV_16U ? (1 << 15) : (1 << 16))) )
 699         sumType = CV_32S;
 700     sumType = CV_MAKETYPE( sumType, cn );
 701
 702     Ptr<BaseRowFilter> rowFilter = getRowSumFilter(srcType, sumType, ksize.width, anchor.x );
 703     Ptr<BaseColumnFilter> columnFilter = getColumnSumFilter(sumType,
 704         dstType, ksize.height, anchor.y, normalize ? 1./(ksize.width*ksize.height) : 1);
 705
 706     return Ptr<FilterEngine>(new FilterEngine(Ptr<BaseFilter>(0), rowFilter, columnFilter,
 707            srcType, dstType, sumType, borderType ));
 708 }
 709
 710
 711 void cv::boxFilter( InputArray _src, OutputArray _dst, int ddepth,
 712                 Size ksize, Point anchor,
 713                 bool normalize, int borderType )
 714 {
 715     Mat src = _src.getMat();
 716     int sdepth = src.depth(), cn = src.channels();
 717     if( ddepth < 0 )
 718         ddepth = sdepth;
 719     _dst.create( src.size(), CV_MAKETYPE(ddepth, cn) );
 720     Mat dst = _dst.getMat();
 721     if( borderType != BORDER_CONSTANT && normalize )
 722     {
 723         if( src.rows == 1 )
 724             ksize.height = 1;
 725         if( src.cols == 1 )
 726             ksize.width = 1;
 727     }
 728 #ifdef HAVE_TEGRA_OPTIMIZATION
 729     if ( tegra::box(src, dst, ksize, anchor, normalize, borderType) )
 730         return;
 731 #endif
 732
 733     Ptr<FilterEngine> f = createBoxFilter( src.type(), dst.type(),
 734                         ksize, anchor, normalize, borderType );
 735     f->apply( src, dst );
 736 }
 737
 738 void cv::blur( InputArray src, OutputArray dst,
 739            Size ksize, Point anchor, int borderType )
 740 {
 741     boxFilter( src, dst, -1, ksize, anchor, true, borderType );
 742 }
 743
 744 /****************************************************************************************\
 745                                      Gaussian Blur
 746 \****************************************************************************************/
 747
 748 cv::Mat cv::getGaussianKernel( int n, double sigma, int ktype )
 749 {
 750     const int SMALL_GAUSSIAN_SIZE = 7;
 751     static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] =
 752     {
 753         {1.f},
 754         {0.25f, 0.5f, 0.25f},
 755         {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f},
 756         {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f}
 757     };
 758
 759     const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ?
 760         small_gaussian_tab[n>>1] : 0;
 761
 762     CV_Assert( ktype == CV_32F || ktype == CV_64F );
 763     Mat kernel(n, 1, ktype);
 764     float* cf = (float*)kernel.data;
 765     double* cd = (double*)kernel.data;
 766
 767     double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8;
 768     double scale2X = -0.5/(sigmaX*sigmaX);
 769     double sum = 0;
 770
 771     int i;
 772     for( i = 0; i < n; i++ )
 773     {
 774         double x = i - (n-1)*0.5;
 775         double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x);
 776         if( ktype == CV_32F )
 777         {
 778             cf[i] = (float)t;
 779             sum += cf[i];
 780         }
 781         else
 782         {
 783             cd[i] = t;
 784             sum += cd[i];
 785         }
 786     }
 787
 788     sum = 1./sum;
 789     for( i = 0; i < n; i++ )
 790     {
 791         if( ktype == CV_32F )
 792             cf[i] = (float)(cf[i]*sum);
 793         else
 794             cd[i] *= sum;
 795     }
 796
 797     return kernel;
 798 }
 799
 800
 801 cv::Ptr<cv::FilterEngine> cv::createGaussianFilter( int type, Size ksize,
 802                                         double sigma1, double sigma2,
 803                                         int borderType )
 804 {
 805     int depth = CV_MAT_DEPTH(type);
 806     if( sigma2 <= 0 )
 807         sigma2 = sigma1;
 808
 809     // automatic detection of kernel size from sigma
 810     if( ksize.width <= 0 && sigma1 > 0 )
 811         ksize.width = cvRound(sigma1*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
 812     if( ksize.height <= 0 && sigma2 > 0 )
 813         ksize.height = cvRound(sigma2*(depth == CV_8U ? 3 : 4)*2 + 1)|1;
 814
 815     CV_Assert( ksize.width > 0 && ksize.width % 2 == 1 &&
 816         ksize.height > 0 && ksize.height % 2 == 1 );
 817
 818     sigma1 = std::max( sigma1, 0. );
 819     sigma2 = std::max( sigma2, 0. );
 820
 821     Mat kx = getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F) );
 822     Mat ky;
 823     if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON )
 824         ky = kx;
 825     else
 826         ky = getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F) );
 827
 828     return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType );
 829 }
 830
 831
 832 void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize,
 833                    double sigma1, double sigma2,
 834                    int borderType )
 835 {
 836     Mat src = _src.getMat();
 837     _dst.create( src.size(), src.type() );
 838     Mat dst = _dst.getMat();
 839
 840     if( borderType != BORDER_CONSTANT )
 841     {
 842         if( src.rows == 1 )
 843             ksize.height = 1;
 844         if( src.cols == 1 )
 845             ksize.width = 1;
 846     }
 847
 848     if( ksize.width == 1 && ksize.height == 1 )
 849     {
 850         src.copyTo(dst);
 851         return;
 852     }
 853
 854 #ifdef HAVE_TEGRA_OPTIMIZATION
 855     if(sigma1 == 0 && sigma2 == 0 && tegra::gaussian(src, dst, ksize, borderType))
 856         return;
 857 #endif
 858
 859     Ptr<FilterEngine> f = createGaussianFilter( src.type(), ksize, sigma1, sigma2, borderType );
 860     f->apply( src, dst );
 861 }
 862
 863
 864 /****************************************************************************************\
 865                                       Median Filter
 866 \****************************************************************************************/
 867
 868 namespace cv
 869 {
 870 typedef ushort HT;
 871
 872 /**
 873  * This structure represents a two-tier histogram. The first tier (known as the
 874  * "coarse" level) is 4 bit wide and the second tier (known as the "fine" level)
 875  * is 8 bit wide. Pixels inserted in the fine level also get inserted into the
 876  * coarse bucket designated by the 4 MSBs of the fine bucket value.
 877  *
 878  * The structure is aligned on 16 bits, which is a prerequisite for SIMD
 879  * instructions. Each bucket is 16 bit wide, which means that extra care must be
 880  * taken to prevent overflow.
 881  */
 882 typedef struct
 883 {
 884     HT coarse[16];
 885     HT fine[16][16];
 886 } Histogram;
 887
 888
 889 #if CV_SSE2
 890 #define MEDIAN_HAVE_SIMD 1
 891
 892 static inline void histogram_add_simd( const HT x[16], HT y[16] )
 893 {
 894     const __m128i* rx = (const __m128i*)x;
 895     __m128i* ry = (__m128i*)y;
 896     __m128i r0 = _mm_add_epi16(_mm_load_si128(ry+0),_mm_load_si128(rx+0));
 897     __m128i r1 = _mm_add_epi16(_mm_load_si128(ry+1),_mm_load_si128(rx+1));
 898     _mm_store_si128(ry+0, r0);
 899     _mm_store_si128(ry+1, r1);
 900 }
 901
 902 static inline void histogram_sub_simd( const HT x[16], HT y[16] )
 903 {
 904     const __m128i* rx = (const __m128i*)x;
 905     __m128i* ry = (__m128i*)y;
 906     __m128i r0 = _mm_sub_epi16(_mm_load_si128(ry+0),_mm_load_si128(rx+0));
 907     __m128i r1 = _mm_sub_epi16(_mm_load_si128(ry+1),_mm_load_si128(rx+1));
 908     _mm_store_si128(ry+0, r0);
 909     _mm_store_si128(ry+1, r1);
 910 }
 911
 912 #else
 913 #define MEDIAN_HAVE_SIMD 0
 914 #endif
 915
 916
 917 static inline void histogram_add( const HT x[16], HT y[16] )
 918 {
 919     int i;
 920     for( i = 0; i < 16; ++i )
 921         y[i] = (HT)(y[i] + x[i]);
 922 }
 923
 924 static inline void histogram_sub( const HT x[16], HT y[16] )
 925 {
 926     int i;
 927     for( i = 0; i < 16; ++i )
 928         y[i] = (HT)(y[i] - x[i]);
 929 }
 930
 931 static inline void histogram_muladd( int a, const HT x[16],
 932         HT y[16] )
 933 {
 934     for( int i = 0; i < 16; ++i )
 935         y[i] = (HT)(y[i] + a * x[i]);
 936 }
 937
 938 static void
 939 medianBlur_8u_O1( const Mat& _src, Mat& _dst, int ksize )
 940 {
 941 /**
 942  * HOP is short for Histogram OPeration. This macro makes an operation \a op on
 943  * histogram \a h for pixel value \a x. It takes care of handling both levels.
 944  */
 945 #define HOP(h,x,op) \
 946     h.coarse[x>>4] op, \
 947     *((HT*)h.fine + x) op
 948
 949 #define COP(c,j,x,op) \
 950     h_coarse[ 16*(n*c+j) + (x>>4) ] op, \
 951     h_fine[ 16 * (n*(16*c+(x>>4)) + j) + (x & 0xF) ] op
 952
 953     int cn = _dst.channels(), m = _dst.rows, r = (ksize-1)/2;
 954     size_t sstep = _src.step, dstep = _dst.step;
 955     Histogram CV_DECL_ALIGNED(16) H[4];
 956     HT CV_DECL_ALIGNED(16) luc[4][16];
 957
 958     int STRIPE_SIZE = std::min( _dst.cols, 512/cn );
 959
 960     vector<HT> _h_coarse(1 * 16 * (STRIPE_SIZE + 2*r) * cn + 16);
 961     vector<HT> _h_fine(16 * 16 * (STRIPE_SIZE + 2*r) * cn + 16);
 962     HT* h_coarse = alignPtr(&_h_coarse[0], 16);
 963     HT* h_fine = alignPtr(&_h_fine[0], 16);
 964 #if MEDIAN_HAVE_SIMD
 965     volatile bool useSIMD = checkHardwareSupport(CV_CPU_SSE2);
 966 #endif
 967
 968     for( int x = 0; x < _dst.cols; x += STRIPE_SIZE )
 969     {
 970         int i, j, k, c, n = std::min(_dst.cols - x, STRIPE_SIZE) + r*2;
 971         const uchar* src = _src.data + x*cn;
 972         uchar* dst = _dst.data + (x - r)*cn;
 973
 974         memset( h_coarse, 0, 16*n*cn*sizeof(h_coarse[0]) );
 975         memset( h_fine, 0, 16*16*n*cn*sizeof(h_fine[0]) );
 976
 977         // First row initialization
 978         for( c = 0; c < cn; c++ )
 979         {
 980             for( j = 0; j < n; j++ )
 981                 COP( c, j, src[cn*j+c], += (cv::HT)(r+2) );
 982
 983             for( i = 1; i < r; i++ )
 984             {
 985                 const uchar* p = src + sstep*std::min(i, m-1);
 986                 for ( j = 0; j < n; j++ )
 987                     COP( c, j, p[cn*j+c], ++ );
 988             }
 989         }
 990
 991         for( i = 0; i < m; i++ )
 992         {
 993             const uchar* p0 = src + sstep * std::max( 0, i-r-1 );
 994             const uchar* p1 = src + sstep * std::min( m-1, i+r );
 995
 996             memset( H, 0, cn*sizeof(H[0]) );
 997             memset( luc, 0, cn*sizeof(luc[0]) );
 998             for( c = 0; c < cn; c++ )
 999             {
1000                 // Update column histograms for the entire row.
1001                 for( j = 0; j < n; j++ )
1002                 {
1003                     COP( c, j, p0[j*cn + c], -- );
1004                     COP( c, j, p1[j*cn + c], ++ );
1005                 }
1006
1007                 // First column initialization
1008                 for( k = 0; k < 16; ++k )
1009                     histogram_muladd( 2*r+1, &h_fine[16*n*(16*c+k)], &H[c].fine[k][0] );
1010
1011             #if MEDIAN_HAVE_SIMD
1012                 if( useSIMD )
1013                 {
1014                     for( j = 0; j < 2*r; ++j )
1015                         histogram_add_simd( &h_coarse[16*(n*c+j)], H[c].coarse );
1016
1017                     for( j = r; j < n-r; j++ )
1018                     {
1019                         int t = 2*r*r + 2*r, b, sum = 0;
1020                         HT* segment;
1021
1022                         histogram_add_simd( &h_coarse[16*(n*c + std::min(j+r,n-1))], H[c].coarse );
1023
1024                         // Find median at coarse level
1025                         for ( k = 0; k < 16 ; ++k )
1026                         {
1027                             sum += H[c].coarse[k];
1028                             if ( sum > t )
1029                             {
1030                                 sum -= H[c].coarse[k];
1031                                 break;
1032                             }
1033                         }
1034                         assert( k < 16 );
1035
1036                         /* Update corresponding histogram segment */
1037                         if ( luc[c][k] <= j-r )
1038                         {
1039                             memset( &H[c].fine[k], 0, 16 * sizeof(HT) );
1040                             for ( luc[c][k] = cv::HT(j-r); luc[c][k] < MIN(j+r+1,n); ++luc[c][k] )
1041                                 histogram_add_simd( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );
1042
1043                             if ( luc[c][k] < j+r+1 )
1044                             {
1045                                 histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );
1046                                 luc[c][k] = (HT)(j+r+1);
1047                             }
1048                         }
1049                         else
1050                         {
1051                             for ( ; luc[c][k] < j+r+1; ++luc[c][k] )
1052                             {
1053                                 histogram_sub_simd( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );
1054                                 histogram_add_simd( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );
1055                             }
1056                         }
1057
1058                         histogram_sub_simd( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );
1059
1060                         /* Find median in segment */
1061                         segment = H[c].fine[k];
1062                         for ( b = 0; b < 16 ; b++ )
1063                         {
1064                             sum += segment[b];
1065                             if ( sum > t )
1066                             {
1067                                 dst[dstep*i+cn*j+c] = (uchar)(16*k + b);
1068                                 break;
1069                             }
1070                         }
1071                         assert( b < 16 );
1072                     }
1073                 }
1074                 else
1075             #endif
1076                 {
1077                     for( j = 0; j < 2*r; ++j )
1078                         histogram_add( &h_coarse[16*(n*c+j)], H[c].coarse );
1079
1080                     for( j = r; j < n-r; j++ )
1081                     {
1082                         int t = 2*r*r + 2*r, b, sum = 0;
1083                         HT* segment;
1084
1085                         histogram_add( &h_coarse[16*(n*c + std::min(j+r,n-1))], H[c].coarse );
1086
1087                         // Find median at coarse level
1088                         for ( k = 0; k < 16 ; ++k )
1089                         {
1090                             sum += H[c].coarse[k];
1091                             if ( sum > t )
1092                             {
1093                                 sum -= H[c].coarse[k];
1094                                 break;
1095                             }
1096                         }
1097                         assert( k < 16 );
1098
1099                         /* Update corresponding histogram segment */
1100                         if ( luc[c][k] <= j-r )
1101                         {
1102                             memset( &H[c].fine[k], 0, 16 * sizeof(HT) );
1103                             for ( luc[c][k] = cv::HT(j-r); luc[c][k] < MIN(j+r+1,n); ++luc[c][k] )
1104                                 histogram_add( &h_fine[16*(n*(16*c+k)+luc[c][k])], H[c].fine[k] );
1105
1106                             if ( luc[c][k] < j+r+1 )
1107                             {
1108                                 histogram_muladd( j+r+1 - n, &h_fine[16*(n*(16*c+k)+(n-1))], &H[c].fine[k][0] );
1109                                 luc[c][k] = (HT)(j+r+1);
1110                             }
1111                         }
1112                         else
1113                         {
1114                             for ( ; luc[c][k] < j+r+1; ++luc[c][k] )
1115                             {
1116                                 histogram_sub( &h_fine[16*(n*(16*c+k)+MAX(luc[c][k]-2*r-1,0))], H[c].fine[k] );
1117                                 histogram_add( &h_fine[16*(n*(16*c+k)+MIN(luc[c][k],n-1))], H[c].fine[k] );
1118                             }
1119                         }
1120
1121                         histogram_sub( &h_coarse[16*(n*c+MAX(j-r,0))], H[c].coarse );
1122
1123                         /* Find median in segment */
1124                         segment = H[c].fine[k];
1125                         for ( b = 0; b < 16 ; b++ )
1126                         {
1127                             sum += segment[b];
1128                             if ( sum > t )
1129                             {
1130                                 dst[dstep*i+cn*j+c] = (uchar)(16*k + b);
1131                                 break;
1132                             }
1133                         }
1134                         assert( b < 16 );
1135                     }
1136                 }
1137             }
1138         }
1139     }
1140
1141 #undef HOP
1142 #undef COP
1143 }
1144
1145 static void
1146 medianBlur_8u_Om( const Mat& _src, Mat& _dst, int m )
1147 {
1148     #define N  16
1149     int     zone0[4][N];
1150     int     zone1[4][N*N];
1151     int     x, y;
1152     int     n2 = m*m/2;
1153     Size    size = _dst.size();
1154     const uchar* src = _src.data;
1155     uchar*  dst = _dst.data;
1156     int     src_step = (int)_src.step, dst_step = (int)_dst.step;
1157     int     cn = _src.channels();
1158     const uchar*  src_max = src + size.height*src_step;
1159
1160     #define UPDATE_ACC01( pix, cn, op ) \
1161     {                                   \
1162         int p = (pix);                  \
1163         zone1[cn][p] op;                \
1164         zone0[cn][p >> 4] op;           \
1165     }
1166
1167     //CV_Assert( size.height >= nx && size.width >= nx );
1168     for( x = 0; x < size.width; x++, src += cn, dst += cn )
1169     {
1170         uchar* dst_cur = dst;
1171         const uchar* src_top = src;
1172         const uchar* src_bottom = src;
1173         int k, c;
1174         int src_step1 = src_step, dst_step1 = dst_step;
1175
1176         if( x % 2 != 0 )
1177         {
1178             src_bottom = src_top += src_step*(size.height-1);
1179             dst_cur += dst_step*(size.height-1);
1180             src_step1 = -src_step1;
1181             dst_step1 = -dst_step1;
1182         }
1183
1184         // init accumulator
1185         memset( zone0, 0, sizeof(zone0[0])*cn );
1186         memset( zone1, 0, sizeof(zone1[0])*cn );
1187
1188         for( y = 0; y <= m/2; y++ )
1189         {
1190             for( c = 0; c < cn; c++ )
1191             {
1192                 if( y > 0 )
1193                 {
1194                     for( k = 0; k < m*cn; k += cn )
1195                         UPDATE_ACC01( src_bottom[k+c], c, ++ );
1196                 }
1197                 else
1198                 {
1199                     for( k = 0; k < m*cn; k += cn )
1200                         UPDATE_ACC01( src_bottom[k+c], c, += m/2+1 );
1201                 }
1202             }
1203
1204             if( (src_step1 > 0 && y < size.height-1) ||
1205                 (src_step1 < 0 && size.height-y-1 > 0) )
1206                 src_bottom += src_step1;
1207         }
1208
1209         for( y = 0; y < size.height; y++, dst_cur += dst_step1 )
1210         {
1211             // find median
1212             for( c = 0; c < cn; c++ )
1213             {
1214                 int s = 0;
1215                 for( k = 0; ; k++ )
1216                 {
1217                     int t = s + zone0[c][k];
1218                     if( t > n2 ) break;
1219                     s = t;
1220                 }
1221
1222                 for( k *= N; ;k++ )
1223                 {
1224                     s += zone1[c][k];
1225                     if( s > n2 ) break;
1226                 }
1227
1228                 dst_cur[c] = (uchar)k;
1229             }
1230
1231             if( y+1 == size.height )
1232                 break;
1233
1234             if( cn == 1 )
1235             {
1236                 for( k = 0; k < m; k++ )
1237                 {
1238                     int p = src_top[k];
1239                     int q = src_bottom[k];
1240                     zone1[0][p]--;
1241                     zone0[0][p>>4]--;
1242                     zone1[0][q]++;
1243                     zone0[0][q>>4]++;
1244                 }
1245             }
1246             else if( cn == 3 )
1247             {
1248                 for( k = 0; k < m*3; k += 3 )
1249                 {
1250                     UPDATE_ACC01( src_top[k], 0, -- );
1251                     UPDATE_ACC01( src_top[k+1], 1, -- );
1252                     UPDATE_ACC01( src_top[k+2], 2, -- );
1253
1254                     UPDATE_ACC01( src_bottom[k], 0, ++ );
1255                     UPDATE_ACC01( src_bottom[k+1], 1, ++ );
1256                     UPDATE_ACC01( src_bottom[k+2], 2, ++ );
1257                 }
1258             }
1259             else
1260             {
1261                 assert( cn == 4 );
1262                 for( k = 0; k < m*4; k += 4 )
1263                 {
1264                     UPDATE_ACC01( src_top[k], 0, -- );
1265                     UPDATE_ACC01( src_top[k+1], 1, -- );
1266                     UPDATE_ACC01( src_top[k+2], 2, -- );
1267                     UPDATE_ACC01( src_top[k+3], 3, -- );
1268
1269                     UPDATE_ACC01( src_bottom[k], 0, ++ );
1270                     UPDATE_ACC01( src_bottom[k+1], 1, ++ );
1271                     UPDATE_ACC01( src_bottom[k+2], 2, ++ );
1272                     UPDATE_ACC01( src_bottom[k+3], 3, ++ );
1273                 }
1274             }
1275
1276             if( (src_step1 > 0 && src_bottom + src_step1 < src_max) ||
1277                 (src_step1 < 0 && src_bottom + src_step1 >= src) )
1278                 src_bottom += src_step1;
1279
1280             if( y >= m/2 )
1281                 src_top += src_step1;
1282         }
1283     }
1284 #undef N
1285 #undef UPDATE_ACC
1286 }
1287
1288
1289 struct MinMax8u
1290 {
1291     typedef uchar value_type;
1292     typedef int arg_type;
1293     enum { SIZE = 1 };
1294     arg_type load(const uchar* ptr) { return *ptr; }
1295     void store(uchar* ptr, arg_type val) { *ptr = (uchar)val; }
1296     void operator()(arg_type& a, arg_type& b) const
1297     {
1298         int t = CV_FAST_CAST_8U(a - b);
1299         b += t; a -= t;
1300     }
1301 };
1302
1303 struct MinMax16u
1304 {
1305     typedef ushort value_type;
1306     typedef int arg_type;
1307     enum { SIZE = 1 };
1308     arg_type load(const ushort* ptr) { return *ptr; }
1309     void store(ushort* ptr, arg_type val) { *ptr = (ushort)val; }
1310     void operator()(arg_type& a, arg_type& b) const
1311     {
1312         arg_type t = a;
1313         a = std::min(a, b);
1314         b = std::max(b, t);
1315     }
1316 };
1317
1318 struct MinMax16s
1319 {
1320     typedef short value_type;
1321     typedef int arg_type;
1322     enum { SIZE = 1 };
1323     arg_type load(const short* ptr) { return *ptr; }
1324     void store(short* ptr, arg_type val) { *ptr = (short)val; }
1325     void operator()(arg_type& a, arg_type& b) const
1326     {
1327         arg_type t = a;
1328         a = std::min(a, b);
1329         b = std::max(b, t);
1330     }
1331 };
1332
1333 struct MinMax32f
1334 {
1335     typedef float value_type;
1336     typedef float arg_type;
1337     enum { SIZE = 1 };
1338     arg_type load(const float* ptr) { return *ptr; }
1339     void store(float* ptr, arg_type val) { *ptr = val; }
1340     void operator()(arg_type& a, arg_type& b) const
1341     {
1342         arg_type t = a;
1343         a = std::min(a, b);
1344         b = std::max(b, t);
1345     }
1346 };
1347
1348 #if CV_SSE2
1349
1350 struct MinMaxVec8u
1351 {
1352     typedef uchar value_type;
1353     typedef __m128i arg_type;
1354     enum { SIZE = 16 };
1355     arg_type load(const uchar* ptr) { return _mm_loadu_si128((const __m128i*)ptr); }
1356     void store(uchar* ptr, arg_type val) { _mm_storeu_si128((__m128i*)ptr, val); }
1357     void operator()(arg_type& a, arg_type& b) const
1358     {
1359         arg_type t = a;
1360         a = _mm_min_epu8(a, b);
1361         b = _mm_max_epu8(b, t);
1362     }
1363 };
1364
1365
1366 struct MinMaxVec16u
1367 {
1368     typedef ushort value_type;
1369     typedef __m128i arg_type;
1370     enum { SIZE = 8 };
1371     arg_type load(const ushort* ptr) { return _mm_loadu_si128((const __m128i*)ptr); }
1372     void store(ushort* ptr, arg_type val) { _mm_storeu_si128((__m128i*)ptr, val); }
1373     void operator()(arg_type& a, arg_type& b) const
1374     {
1375         arg_type t = _mm_subs_epu16(a, b);
1376         a = _mm_subs_epu16(a, t);
1377         b = _mm_adds_epu16(b, t);
1378     }
1379 };
1380
1381
1382 struct MinMaxVec16s
1383 {
1384     typedef short value_type;
1385     typedef __m128i arg_type;
1386     enum { SIZE = 8 };
1387     arg_type load(const short* ptr) { return _mm_loadu_si128((const __m128i*)ptr); }
1388     void store(short* ptr, arg_type val) { _mm_storeu_si128((__m128i*)ptr, val); }
1389     void operator()(arg_type& a, arg_type& b) const
1390     {
1391         arg_type t = a;
1392         a = _mm_min_epi16(a, b);
1393         b = _mm_max_epi16(b, t);
1394     }
1395 };
1396
1397
1398 struct MinMaxVec32f
1399 {
1400     typedef float value_type;
1401     typedef __m128 arg_type;
1402     enum { SIZE = 4 };
1403     arg_type load(const float* ptr) { return _mm_loadu_ps(ptr); }
1404     void store(float* ptr, arg_type val) { _mm_storeu_ps(ptr, val); }
1405     void operator()(arg_type& a, arg_type& b) const
1406     {
1407         arg_type t = a;
1408         a = _mm_min_ps(a, b);
1409         b = _mm_max_ps(b, t);
1410     }
1411 };
1412
1413
1414 #else
1415
1416 typedef MinMax8u MinMaxVec8u;
1417 typedef MinMax16u MinMaxVec16u;
1418 typedef MinMax16s MinMaxVec16s;
1419 typedef MinMax32f MinMaxVec32f;
1420
1421 #endif
1422
1423 template<class Op, class VecOp>
1424 static void
1425 medianBlur_SortNet( const Mat& _src, Mat& _dst, int m )
1426 {
1427     typedef typename Op::value_type T;
1428     typedef typename Op::arg_type WT;
1429     typedef typename VecOp::arg_type VT;
1430
1431     const T* src = (const T*)_src.data;
1432     T* dst = (T*)_dst.data;
1433     int sstep = (int)(_src.step/sizeof(T));
1434     int dstep = (int)(_dst.step/sizeof(T));
1435     Size size = _dst.size();
1436     int i, j, k, cn = _src.channels();
1437     Op op;
1438     VecOp vop;
1439     volatile bool useSIMD = checkHardwareSupport(CV_CPU_SSE2);
1440
1441     if( m == 3 )
1442     {
1443         if( size.width == 1 || size.height == 1 )
1444         {
1445             int len = size.width + size.height - 1;
1446             int sdelta = size.height == 1 ? cn : sstep;
1447             int sdelta0 = size.height == 1 ? 0 : sstep - cn;
1448             int ddelta = size.height == 1 ? cn : dstep;
1449
1450             for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
1451                 for( j = 0; j < cn; j++, src++ )
1452                 {
1453                     WT p0 = src[i > 0 ? -sdelta : 0];
1454                     WT p1 = src[0];
1455                     WT p2 = src[i < len - 1 ? sdelta : 0];
1456
1457                     op(p0, p1); op(p1, p2); op(p0, p1);
1458                     dst[j] = (T)p1;
1459                 }
1460             return;
1461         }
1462
1463         size.width *= cn;
1464         for( i = 0; i < size.height; i++, dst += dstep )
1465         {
1466             const T* row0 = src + std::max(i - 1, 0)*sstep;
1467             const T* row1 = src + i*sstep;
1468             const T* row2 = src + std::min(i + 1, size.height-1)*sstep;
1469             int limit = useSIMD ? cn : size.width;
1470
1471             for(j = 0;; )
1472             {
1473                 for( ; j < limit; j++ )
1474                 {
1475                     int j0 = j >= cn ? j - cn : j;
1476                     int j2 = j < size.width - cn ? j + cn : j;
1477                     WT p0 = row0[j0], p1 = row0[j], p2 = row0[j2];
1478                     WT p3 = row1[j0], p4 = row1[j], p5 = row1[j2];
1479                     WT p6 = row2[j0], p7 = row2[j], p8 = row2[j2];
1480
1481                     op(p1, p2); op(p4, p5); op(p7, p8); op(p0, p1);
1482                     op(p3, p4); op(p6, p7); op(p1, p2); op(p4, p5);
1483                     op(p7, p8); op(p0, p3); op(p5, p8); op(p4, p7);
1484                     op(p3, p6); op(p1, p4); op(p2, p5); op(p4, p7);
1485                     op(p4, p2); op(p6, p4); op(p4, p2);
1486                     dst[j] = (T)p4;
1487                 }
1488
1489                 if( limit == size.width )
1490                     break;
1491
1492                 for( ; j <= size.width - VecOp::SIZE - cn; j += VecOp::SIZE )
1493                 {
1494                     VT p0 = vop.load(row0+j-cn), p1 = vop.load(row0+j), p2 = vop.load(row0+j+cn);
1495                     VT p3 = vop.load(row1+j-cn), p4 = vop.load(row1+j), p5 = vop.load(row1+j+cn);
1496                     VT p6 = vop.load(row2+j-cn), p7 = vop.load(row2+j), p8 = vop.load(row2+j+cn);
1497
1498                     vop(p1, p2); vop(p4, p5); vop(p7, p8); vop(p0, p1);
1499                     vop(p3, p4); vop(p6, p7); vop(p1, p2); vop(p4, p5);
1500                     vop(p7, p8); vop(p0, p3); vop(p5, p8); vop(p4, p7);
1501                     vop(p3, p6); vop(p1, p4); vop(p2, p5); vop(p4, p7);
1502                     vop(p4, p2); vop(p6, p4); vop(p4, p2);
1503                     vop.store(dst+j, p4);
1504                 }
1505
1506                 limit = size.width;
1507             }
1508         }
1509     }
1510     else if( m == 5 )
1511     {
1512         if( size.width == 1 || size.height == 1 )
1513         {
1514             int len = size.width + size.height - 1;
1515             int sdelta = size.height == 1 ? cn : sstep;
1516             int sdelta0 = size.height == 1 ? 0 : sstep - cn;
1517             int ddelta = size.height == 1 ? cn : dstep;
1518
1519             for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
1520                 for( j = 0; j < cn; j++, src++ )
1521                 {
1522                     int i1 = i > 0 ? -sdelta : 0;
1523                     int i0 = i > 1 ? -sdelta*2 : i1;
1524                     int i3 = i < len-1 ? sdelta : 0;
1525                     int i4 = i < len-2 ? sdelta*2 : i3;
1526                     WT p0 = src[i0], p1 = src[i1], p2 = src[0], p3 = src[i3], p4 = src[i4];
1527
1528                     op(p0, p1); op(p3, p4); op(p2, p3); op(p3, p4); op(p0, p2);
1529                     op(p2, p4); op(p1, p3); op(p1, p2);
1530                     dst[j] = (T)p2;
1531                 }
1532             return;
1533         }
1534
1535         size.width *= cn;
1536         for( i = 0; i < size.height; i++, dst += dstep )
1537         {
1538             const T* row[5];
1539             row[0] = src + std::max(i - 2, 0)*sstep;
1540             row[1] = src + std::max(i - 1, 0)*sstep;
1541             row[2] = src + i*sstep;
1542             row[3] = src + std::min(i + 1, size.height-1)*sstep;
1543             row[4] = src + std::min(i + 2, size.height-1)*sstep;
1544             int limit = useSIMD ? cn*2 : size.width;
1545
1546             for(j = 0;; )
1547             {
1548                 for( ; j < limit; j++ )
1549                 {
1550                     WT p[25];
1551                     int j1 = j >= cn ? j - cn : j;
1552                     int j0 = j >= cn*2 ? j - cn*2 : j1;
1553                     int j3 = j < size.width - cn ? j + cn : j;
1554                     int j4 = j < size.width - cn*2 ? j + cn*2 : j3;
1555                     for( k = 0; k < 5; k++ )
1556                     {
1557                         const T* rowk = row[k];
1558                         p[k*5] = rowk[j0]; p[k*5+1] = rowk[j1];
1559                         p[k*5+2] = rowk[j]; p[k*5+3] = rowk[j3];
1560                         p[k*5+4] = rowk[j4];
1561                     }
1562
1563                     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]);
1564                     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]);
1565                     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]);
1566                     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]);
1567                     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]);
1568                     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]);
1569                     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]);
1570                     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]);
1571                     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]);
1572                     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]);
1573                     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]);
1574                     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]);
1575                     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]);
1576                     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]);
1577                     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]);
1578                     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]);
1579                     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]);
1580                     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]);
1581                     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]);
1582                     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]);
1583                     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]);
1584                     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]);
1585                     op(p[7], p[11]); op(p[11], p[13]); op(p[11], p[12]);
1586                     dst[j] = (T)p[12];
1587                 }
1588
1589                 if( limit == size.width )
1590                     break;
1591
1592                 for( ; j <= size.width - VecOp::SIZE - cn*2; j += VecOp::SIZE )
1593                 {
1594                     VT p[25];
1595                     for( k = 0; k < 5; k++ )
1596                     {
1597                         const T* rowk = row[k];
1598                         p[k*5] = vop.load(rowk+j-cn*2); p[k*5+1] = vop.load(rowk+j-cn);
1599                         p[k*5+2] = vop.load(rowk+j); p[k*5+3] = vop.load(rowk+j+cn);
1600                         p[k*5+4] = vop.load(rowk+j+cn*2);
1601                     }
1602
1603                     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]);
1604                     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]);
1605                     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]);
1606                     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]);
1607                     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]);
1608                     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]);
1609                     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]);
1610                     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]);
1611                     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]);
1612                     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]);
1613                     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]);
1614                     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]);
1615                     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]);
1616                     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]);
1617                     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]);
1618                     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]);
1619                     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]);
1620                     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]);
1621                     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]);
1622                     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]);
1623                     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]);
1624                     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]);
1625                     vop(p[7], p[11]); vop(p[11], p[13]); vop(p[11], p[12]);
1626                     vop.store(dst+j, p[12]);
1627                 }
1628
1629                 limit = size.width;
1630             }
1631         }
1632     }
1633 }
1634
1635 }
1636
1637 void cv::medianBlur( InputArray _src0, OutputArray _dst, int ksize )
1638 {
1639     Mat src0 = _src0.getMat();
1640     _dst.create( src0.size(), src0.type() );
1641     Mat dst = _dst.getMat();
1642
1643     if( ksize <= 1 )
1644     {
1645         src0.copyTo(dst);
1646         return;
1647     }
1648
1649     CV_Assert( ksize % 2 == 1 );
1650
1651 #ifdef HAVE_TEGRA_OPTIMIZATION
1652     if (tegra::medianBlur(src0, dst, ksize))
1653         return;
1654 #endif
1655
1656     bool useSortNet = ksize == 3 || (ksize == 5
1657 #if !CV_SSE2
1658             && src0.depth() > CV_8U
1659 #endif
1660         );
1661
1662     Mat src;
1663     if( useSortNet )
1664     {
1665         if( dst.data != src0.data )
1666             src = src0;
1667         else
1668             src0.copyTo(src);
1669
1670         if( src.depth() == CV_8U )
1671             medianBlur_SortNet<MinMax8u, MinMaxVec8u>( src, dst, ksize );
1672         else if( src.depth() == CV_16U )
1673             medianBlur_SortNet<MinMax16u, MinMaxVec16u>( src, dst, ksize );
1674         else if( src.depth() == CV_16S )
1675             medianBlur_SortNet<MinMax16s, MinMaxVec16s>( src, dst, ksize );
1676         else if( src.depth() == CV_32F )
1677             medianBlur_SortNet<MinMax32f, MinMaxVec32f>( src, dst, ksize );
1678         else
1679             CV_Error(CV_StsUnsupportedFormat, "");
1680
1681         return;
1682     }
1683     else
1684     {
1685         cv::copyMakeBorder( src0, src, 0, 0, ksize/2, ksize/2, BORDER_REPLICATE );
1686
1687         int cn = src0.channels();
1688         CV_Assert( src.depth() == CV_8U && (cn == 1 || cn == 3 || cn == 4) );
1689
1690         double img_size_mp = (double)(src0.total())/(1 << 20);
1691         if( ksize <= 3 + (img_size_mp < 1 ? 12 : img_size_mp < 4 ? 6 : 2)*(MEDIAN_HAVE_SIMD && checkHardwareSupport(CV_CPU_SSE2) ? 1 : 3))
1692             medianBlur_8u_Om( src, dst, ksize );
1693         else
1694             medianBlur_8u_O1( src, dst, ksize );
1695     }
1696 }
1697
1698 /****************************************************************************************\
1699                                    Bilateral Filtering
1700 \****************************************************************************************/
1701
1702 namespace cv
1703 {
1704
1705 class BilateralFilter_8u_Invoker :
1706     public ParallelLoopBody
1707 {
1708 public:
1709     BilateralFilter_8u_Invoker(Mat& _dest, const Mat& _temp, int _radius, int _maxk,
1710         int* _space_ofs, float *_space_weight, float *_color_weight) :
1711         temp(&_temp), dest(&_dest), radius(_radius),
1712         maxk(_maxk), space_ofs(_space_ofs), space_weight(_space_weight), color_weight(_color_weight)
1713     {
1714     }
1715
1716     virtual void operator() (const Range& range) const
1717     {
1718         int i, j, cn = dest->channels(), k;
1719         Size size = dest->size();
1720         #if CV_SSE3
1721         int CV_DECL_ALIGNED(16) buf[4];
1722         float CV_DECL_ALIGNED(16) bufSum[4];
1723         static const int CV_DECL_ALIGNED(16) bufSignMask[] = { 0x80000000, 0x80000000, 0x80000000, 0x80000000 };
1724         bool haveSSE3 = checkHardwareSupport(CV_CPU_SSE3);
1725         #endif
1726
1727         for( i = range.start; i < range.end; i++ )
1728         {
1729             const uchar* sptr = temp->ptr(i+radius) + radius*cn;
1730             uchar* dptr = dest->ptr(i);
1731
1732             if( cn == 1 )
1733             {
1734                 for( j = 0; j < size.width; j++ )
1735                 {
1736                     float sum = 0, wsum = 0;
1737                     int val0 = sptr[j];
1738                     k = 0;
1739                     #if CV_SSE3
1740                     if( haveSSE3 )
1741                     {
1742                         __m128 _val0 = _mm_set1_ps(static_cast<float>(val0));
1743                         const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
1744
1745                         for( ; k <= maxk - 4; k += 4 )
1746                         {
1747                             __m128 _valF = _mm_set_ps(sptr[j + space_ofs[k+3]], sptr[j + space_ofs[k+2]],
1748                                                       sptr[j + space_ofs[k+1]], sptr[j + space_ofs[k]]);
1749
1750                             __m128 _val = _mm_andnot_ps(_signMask, _mm_sub_ps(_valF, _val0));
1751                             _mm_store_si128((__m128i*)buf, _mm_cvtps_epi32(_val));
1752
1753                             __m128 _cw = _mm_set_ps(color_weight[buf[3]],color_weight[buf[2]],
1754                                                     color_weight[buf[1]],color_weight[buf[0]]);
1755                             __m128 _sw = _mm_loadu_ps(space_weight+k);
1756                             __m128 _w = _mm_mul_ps(_cw, _sw);
1757                              _cw = _mm_mul_ps(_w, _valF);
1758
1759                              _sw = _mm_hadd_ps(_w, _cw);
1760                              _sw = _mm_hadd_ps(_sw, _sw);
1761                              _mm_storel_pi((__m64*)bufSum, _sw);
1762
1763                              sum += bufSum[1];
1764                              wsum += bufSum[0];
1765                         }
1766                     }
1767                     #endif
1768                     for( ; k < maxk; k++ )
1769                     {
1770                         int val = sptr[j + space_ofs[k]];
1771                         float w = space_weight[k]*color_weight[std::abs(val - val0)];
1772                         sum += val*w;
1773                         wsum += w;
1774                     }
1775                     // overflow is not possible here => there is no need to use CV_CAST_8U
1776                     dptr[j] = (uchar)cvRound(sum/wsum);
1777                 }
1778             }
1779             else
1780             {
1781                 assert( cn == 3 );
1782                 for( j = 0; j < size.width*3; j += 3 )
1783                 {
1784                     float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
1785                     int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
1786                     k = 0;
1787                     #if CV_SSE3
1788                     if( haveSSE3 )
1789                     {
1790                         const __m128 _b0 = _mm_set1_ps(static_cast<float>(b0));
1791                         const __m128 _g0 = _mm_set1_ps(static_cast<float>(g0));
1792                         const __m128 _r0 = _mm_set1_ps(static_cast<float>(r0));
1793                         const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
1794
1795                         for( ; k <= maxk - 4; k += 4 )
1796                         {
1797                             const uchar* sptr_k  = sptr + j + space_ofs[k];
1798                             const uchar* sptr_k1 = sptr + j + space_ofs[k+1];
1799                             const uchar* sptr_k2 = sptr + j + space_ofs[k+2];
1800                             const uchar* sptr_k3 = sptr + j + space_ofs[k+3];
1801
1802                             __m128 _b = _mm_set_ps(sptr_k3[0],sptr_k2[0],sptr_k1[0],sptr_k[0]);
1803                             __m128 _g = _mm_set_ps(sptr_k3[1],sptr_k2[1],sptr_k1[1],sptr_k[1]);
1804                             __m128 _r = _mm_set_ps(sptr_k3[2],sptr_k2[2],sptr_k1[2],sptr_k[2]);
1805
1806                             __m128 bt = _mm_andnot_ps(_signMask, _mm_sub_ps(_b,_b0));
1807                             __m128 gt = _mm_andnot_ps(_signMask, _mm_sub_ps(_g,_g0));
1808                             __m128 rt = _mm_andnot_ps(_signMask, _mm_sub_ps(_r,_r0));
1809
1810                             bt =_mm_add_ps(rt, _mm_add_ps(bt, gt));
1811                             _mm_store_si128((__m128i*)buf, _mm_cvtps_epi32(bt));
1812
1813                             __m128 _w  = _mm_set_ps(color_weight[buf[3]],color_weight[buf[2]],
1814                                                     color_weight[buf[1]],color_weight[buf[0]]);
1815                             __m128 _sw = _mm_loadu_ps(space_weight+k);
1816
1817                             _w = _mm_mul_ps(_w,_sw);
1818                             _b = _mm_mul_ps(_b, _w);
1819                             _g = _mm_mul_ps(_g, _w);
1820                             _r = _mm_mul_ps(_r, _w);
1821
1822                              _w = _mm_hadd_ps(_w, _b);
1823                              _g = _mm_hadd_ps(_g, _r);
1824
1825                              _w = _mm_hadd_ps(_w, _g);
1826                              _mm_store_ps(bufSum, _w);
1827
1828                              wsum  += bufSum[0];
1829                              sum_b += bufSum[1];
1830                              sum_g += bufSum[2];
1831                              sum_r += bufSum[3];
1832                          }
1833                     }
1834                     #endif
1835
1836                     for( ; k < maxk; k++ )
1837                     {
1838                         const uchar* sptr_k = sptr + j + space_ofs[k];
1839                         int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
1840                         float w = space_weight[k]*color_weight[std::abs(b - b0) +
1841                                                                std::abs(g - g0) + std::abs(r - r0)];
1842                         sum_b += b*w; sum_g += g*w; sum_r += r*w;
1843                         wsum += w;
1844                     }
1845                     wsum = 1.f/wsum;
1846                     b0 = cvRound(sum_b*wsum);
1847                     g0 = cvRound(sum_g*wsum);
1848                     r0 = cvRound(sum_r*wsum);
1849                     dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0;
1850                 }
1851             }
1852         }
1853     }
1854
1855 private:
1856     const Mat *temp;
1857     Mat *dest;
1858     int radius, maxk, *space_ofs;
1859     float *space_weight, *color_weight;
1860 };
1861
1862 static void
1863 bilateralFilter_8u( const Mat& src, Mat& dst, int d,
1864     double sigma_color, double sigma_space,
1865     int borderType )
1866 {
1867
1868     int cn = src.channels();
1869     int i, j, maxk, radius;
1870     Size size = src.size();
1871
1872     CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) &&
1873               src.type() == dst.type() && src.size() == dst.size() &&
1874               src.data != dst.data );
1875
1876     if( sigma_color <= 0 )
1877         sigma_color = 1;
1878     if( sigma_space <= 0 )
1879         sigma_space = 1;
1880
1881     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
1882     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
1883
1884     if( d <= 0 )
1885         radius = cvRound(sigma_space*1.5);
1886     else
1887         radius = d/2;
1888     radius = MAX(radius, 1);
1889     d = radius*2 + 1;
1890
1891     Mat temp;
1892     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
1893
1894     vector<float> _color_weight(cn*256);
1895     vector<float> _space_weight(d*d);
1896     vector<int> _space_ofs(d*d);
1897     float* color_weight = &_color_weight[0];
1898     float* space_weight = &_space_weight[0];
1899     int* space_ofs = &_space_ofs[0];
1900
1901     // initialize color-related bilateral filter coefficients
1902
1903     for( i = 0; i < 256*cn; i++ )
1904         color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);
1905
1906     // initialize space-related bilateral filter coefficients
1907     for( i = -radius, maxk = 0; i <= radius; i++ )
1908     {
1909         j = -radius;
1910
1911         for( ;j <= radius; j++ )
1912         {
1913             double r = std::sqrt((double)i*i + (double)j*j);
1914             if( r > radius )
1915                 continue;
1916             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
1917             space_ofs[maxk++] = (int)(i*temp.step + j*cn);
1918         }
1919     }
1920
1921     BilateralFilter_8u_Invoker body(dst, temp, radius, maxk, space_ofs, space_weight, color_weight);
1922     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
1923 }
1924
1925
1926 class BilateralFilter_32f_Invoker :
1927     public ParallelLoopBody
1928 {
1929 public:
1930
1931     BilateralFilter_32f_Invoker(int _cn, int _radius, int _maxk, int *_space_ofs,
1932         const Mat& _temp, Mat& _dest, float _scale_index, float *_space_weight, float *_expLUT) :
1933         cn(_cn), radius(_radius), maxk(_maxk), space_ofs(_space_ofs),
1934         temp(&_temp), dest(&_dest), scale_index(_scale_index), space_weight(_space_weight), expLUT(_expLUT)
1935     {
1936     }
1937
1938     virtual void operator() (const Range& range) const
1939     {
1940         int i, j, k;
1941         Size size = dest->size();
1942         #if CV_SSE3
1943         int CV_DECL_ALIGNED(16) idxBuf[4];
1944         float CV_DECL_ALIGNED(16) bufSum32[4];
1945         static const int CV_DECL_ALIGNED(16) bufSignMask[] = { 0x80000000, 0x80000000, 0x80000000, 0x80000000 };
1946         bool haveSSE3 = checkHardwareSupport(CV_CPU_SSE3);
1947         #endif
1948
1949         for( i = range.start; i < range.end; i++ )
1950         {
1951             const float* sptr = temp->ptr<float>(i+radius) + radius*cn;
1952             float* dptr = dest->ptr<float>(i);
1953
1954             if( cn == 1 )
1955             {
1956                 for( j = 0; j < size.width; j++ )
1957                 {
1958                     float sum = 0, wsum = 0;
1959                     float val0 = sptr[j];
1960                     k = 0;
1961                     #if CV_SSE3
1962                     if( haveSSE3 )
1963                     {
1964                         const __m128 _val0 = _mm_set1_ps(sptr[j]);
1965                         const __m128 _scale_index = _mm_set1_ps(scale_index);
1966                         const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
1967
1968                         for( ; k <= maxk - 4 ; k += 4 )
1969                         {
1970                             __m128 _sw    = _mm_loadu_ps(space_weight + k);
1971                             __m128 _val   = _mm_set_ps(sptr[j + space_ofs[k+3]], sptr[j + space_ofs[k+2]],
1972                                                        sptr[j + space_ofs[k+1]], sptr[j + space_ofs[k]]);
1973                             __m128 _alpha = _mm_mul_ps(_mm_andnot_ps( _signMask, _mm_sub_ps(_val,_val0)), _scale_index);
1974
1975                             __m128i _idx = _mm_cvtps_epi32(_alpha);
1976                             _mm_store_si128((__m128i*)idxBuf, _idx);
1977                             _alpha = _mm_sub_ps(_alpha, _mm_cvtepi32_ps(_idx));
1978
1979                             __m128 _explut  = _mm_set_ps(expLUT[idxBuf[3]], expLUT[idxBuf[2]],
1980                                                          expLUT[idxBuf[1]], expLUT[idxBuf[0]]);
1981                             __m128 _explut1 = _mm_set_ps(expLUT[idxBuf[3]+1], expLUT[idxBuf[2]+1],
1982                                                          expLUT[idxBuf[1]+1], expLUT[idxBuf[0]+1]);
1983
1984                             __m128 _w = _mm_mul_ps(_sw, _mm_add_ps(_explut, _mm_mul_ps(_alpha, _mm_sub_ps(_explut1, _explut))));
1985                             _val = _mm_mul_ps(_w, _val);
1986
1987                              _sw = _mm_hadd_ps(_w, _val);
1988                              _sw = _mm_hadd_ps(_sw, _sw);
1989                              _mm_storel_pi((__m64*)bufSum32, _sw);
1990
1991                              sum += bufSum32[1];
1992                              wsum += bufSum32[0];
1993                         }
1994                     }
1995                     #endif
1996
1997                     for( ; k < maxk; k++ )
1998                     {
1999                         float val = sptr[j + space_ofs[k]];
2000                         float alpha = (float)(std::abs(val - val0)*scale_index);
2001                         int idx = cvFloor(alpha);
2002                         alpha -= idx;
2003                         float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
2004                         sum += val*w;
2005                         wsum += w;
2006                     }
2007                     dptr[j] = (float)(sum/wsum);
2008                 }
2009             }
2010             else
2011             {
2012                 assert( cn == 3 );
2013                 for( j = 0; j < size.width*3; j += 3 )
2014                 {
2015                     float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
2016                     float b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
2017                     k = 0;
2018                     #if  CV_SSE3
2019                     if( haveSSE3 )
2020                     {
2021                         const __m128 _b0 = _mm_set1_ps(b0);
2022                         const __m128 _g0 = _mm_set1_ps(g0);
2023                         const __m128 _r0 = _mm_set1_ps(r0);
2024                         const __m128 _scale_index = _mm_set1_ps(scale_index);
2025                         const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
2026
2027                         for( ; k <= maxk-4; k += 4 )
2028                         {
2029                             __m128 _sw = _mm_loadu_ps(space_weight + k);
2030
2031                             const float* sptr_k  = sptr + j + space_ofs[k];
2032                             const float* sptr_k1 = sptr + j + space_ofs[k+1];
2033                             const float* sptr_k2 = sptr + j + space_ofs[k+2];
2034                             const float* sptr_k3 = sptr + j + space_ofs[k+3];
2035
2036                             __m128 _b = _mm_set_ps(sptr_k3[0], sptr_k2[0], sptr_k1[0], sptr_k[0]);
2037                             __m128 _g = _mm_set_ps(sptr_k3[1], sptr_k2[1], sptr_k1[1], sptr_k[1]);
2038                             __m128 _r = _mm_set_ps(sptr_k3[2], sptr_k2[2], sptr_k1[2], sptr_k[2]);
2039
2040                             __m128 _bt = _mm_andnot_ps(_signMask,_mm_sub_ps(_b,_b0));
2041                             __m128 _gt = _mm_andnot_ps(_signMask,_mm_sub_ps(_g,_g0));
2042                             __m128 _rt = _mm_andnot_ps(_signMask,_mm_sub_ps(_r,_r0));
2043
2044                             __m128 _alpha = _mm_mul_ps(_scale_index, _mm_add_ps(_rt,_mm_add_ps(_bt, _gt)));
2045
2046                             __m128i _idx  = _mm_cvtps_epi32(_alpha);
2047                             _mm_store_si128((__m128i*)idxBuf, _idx);
2048                             _alpha = _mm_sub_ps(_alpha, _mm_cvtepi32_ps(_idx));
2049
2050                             __m128 _explut  = _mm_set_ps(expLUT[idxBuf[3]], expLUT[idxBuf[2]], expLUT[idxBuf[1]], expLUT[idxBuf[0]]);
2051                             __m128 _explut1 = _mm_set_ps(expLUT[idxBuf[3]+1], expLUT[idxBuf[2]+1], expLUT[idxBuf[1]+1], expLUT[idxBuf[0]+1]);
2052
2053                             __m128 _w = _mm_mul_ps(_sw, _mm_add_ps(_explut, _mm_mul_ps(_alpha, _mm_sub_ps(_explut1, _explut))));
2054
2055                             _b = _mm_mul_ps(_b, _w);
2056                             _g = _mm_mul_ps(_g, _w);
2057                             _r = _mm_mul_ps(_r, _w);
2058
2059                              _w = _mm_hadd_ps(_w, _b);
2060                              _g = _mm_hadd_ps(_g, _r);
2061
2062                              _w = _mm_hadd_ps(_w, _g);
2063                              _mm_store_ps(bufSum32, _w);
2064
2065                              wsum  += bufSum32[0];
2066                              sum_b += bufSum32[1];
2067                              sum_g += bufSum32[2];
2068                              sum_r += bufSum32[3];
2069                         }
2070
2071                     }
2072                     #endif
2073
2074                     for(; k < maxk; k++ )
2075                     {
2076                         const float* sptr_k = sptr + j + space_ofs[k];
2077                         float b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
2078                         float alpha = (float)((std::abs(b - b0) +
2079                             std::abs(g - g0) + std::abs(r - r0))*scale_index);
2080                         int idx = cvFloor(alpha);
2081                         alpha -= idx;
2082                         float w = space_weight[k]*(expLUT[idx] + alpha*(expLUT[idx+1] - expLUT[idx]));
2083                         sum_b += b*w; sum_g += g*w; sum_r += r*w;
2084                         wsum += w;
2085                     }
2086                     wsum = 1.f/wsum;
2087                     b0 = sum_b*wsum;
2088                     g0 = sum_g*wsum;
2089                     r0 = sum_r*wsum;
2090                     dptr[j] = b0; dptr[j+1] = g0; dptr[j+2] = r0;
2091                 }
2092             }
2093         }
2094     }
2095
2096 private:
2097     int cn, radius, maxk, *space_ofs;
2098     const Mat* temp;
2099     Mat *dest;
2100     float scale_index, *space_weight, *expLUT;
2101 };
2102
2103
2104 static void
2105 bilateralFilter_32f( const Mat& src, Mat& dst, int d,
2106                      double sigma_color, double sigma_space,
2107                      int borderType )
2108 {
2109     int cn = src.channels();
2110     int i, j, maxk, radius;
2111     double minValSrc=-1, maxValSrc=1;
2112     const int kExpNumBinsPerChannel = 1 << 12;
2113     int kExpNumBins = 0;
2114     float lastExpVal = 1.f;
2115     float len, scale_index;
2116     Size size = src.size();
2117
2118     CV_Assert( (src.type() == CV_32FC1 || src.type() == CV_32FC3) &&
2119         src.type() == dst.type() && src.size() == dst.size() &&
2120         src.data != dst.data );
2121
2122     if( sigma_color <= 0 )
2123         sigma_color = 1;
2124     if( sigma_space <= 0 )
2125         sigma_space = 1;
2126
2127     double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
2128     double gauss_space_coeff = -0.5/(sigma_space*sigma_space);
2129
2130     if( d <= 0 )
2131         radius = cvRound(sigma_space*1.5);
2132     else
2133         radius = d/2;
2134     radius = MAX(radius, 1);
2135     d = radius*2 + 1;
2136     // compute the min/max range for the input image (even if multichannel)
2137
2138     minMaxLoc( src.reshape(1), &minValSrc, &maxValSrc );
2139     if(std::abs(minValSrc - maxValSrc) < FLT_EPSILON)
2140     {
2141         src.copyTo(dst);
2142         return;
2143     }
2144
2145     // temporary copy of the image with borders for easy processing
2146     Mat temp;
2147     copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );
2148     const double insteadNaNValue = -5. * sigma_color;
2149     patchNaNs( temp, insteadNaNValue ); // this replacement of NaNs makes the assumption that depth values are nonnegative
2150                                         // TODO: make insteadNaNValue avalible in the outside function interface to control the cases breaking the assumption
2151     // allocate lookup tables
2152     vector<float> _space_weight(d*d);
2153     vector<int> _space_ofs(d*d);
2154     float* space_weight = &_space_weight[0];
2155     int* space_ofs = &_space_ofs[0];
2156
2157     // assign a length which is slightly more than needed
2158     len = (float)(maxValSrc - minValSrc) * cn;
2159     kExpNumBins = kExpNumBinsPerChannel * cn;
2160     vector<float> _expLUT(kExpNumBins+2);
2161     float* expLUT = &_expLUT[0];
2162
2163     scale_index = kExpNumBins/len;
2164
2165     // initialize the exp LUT
2166     for( i = 0; i < kExpNumBins+2; i++ )
2167     {
2168         if( lastExpVal > 0.f )
2169         {
2170             double val =  i / scale_index;
2171             expLUT[i] = (float)std::exp(val * val * gauss_color_coeff);
2172             lastExpVal = expLUT[i];
2173         }
2174         else
2175             expLUT[i] = 0.f;
2176     }
2177
2178     // initialize space-related bilateral filter coefficients
2179     for( i = -radius, maxk = 0; i <= radius; i++ )
2180         for( j = -radius; j <= radius; j++ )
2181         {
2182             double r = std::sqrt((double)i*i + (double)j*j);
2183             if( r > radius )
2184                 continue;
2185             space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
2186             space_ofs[maxk++] = (int)(i*(temp.step/sizeof(float)) + j*cn);
2187         }
2188
2189     // parallel_for usage
2190
2191     BilateralFilter_32f_Invoker body(cn, radius, maxk, space_ofs, temp, dst, scale_index, space_weight, expLUT);
2192     parallel_for_(Range(0, size.height), body, dst.total()/(double)(1<<16));
2193 }
2194
2195 }
2196
2197 void cv::bilateralFilter( InputArray _src, OutputArray _dst, int d,
2198                       double sigmaColor, double sigmaSpace,
2199                       int borderType )
2200 {
2201     Mat src = _src.getMat();
2202     _dst.create( src.size(), src.type() );
2203     Mat dst = _dst.getMat();
2204
2205     if( src.depth() == CV_8U )
2206         bilateralFilter_8u( src, dst, d, sigmaColor, sigmaSpace, borderType );
2207     else if( src.depth() == CV_32F )
2208         bilateralFilter_32f( src, dst, d, sigmaColor, sigmaSpace, borderType );
2209     else
2210         CV_Error( CV_StsUnsupportedFormat,
2211         "Bilateral filtering is only implemented for 8u and 32f images" );
2212 }
2213
2214 //////////////////////////////////////////////////////////////////////////////////////////
2215
2216 CV_IMPL void
2217 cvSmooth( const void* srcarr, void* dstarr, int smooth_type,
2218           int param1, int param2, double param3, double param4 )
2219 {
2220     cv::Mat src = cv::cvarrToMat(srcarr), dst0 = cv::cvarrToMat(dstarr), dst = dst0;
2221
2222     CV_Assert( dst.size() == src.size() &&
2223         (smooth_type == CV_BLUR_NO_SCALE || dst.type() == src.type()) );
2224
2225     if( param2 <= 0 )
2226         param2 = param1;
2227
2228     if( smooth_type == CV_BLUR || smooth_type == CV_BLUR_NO_SCALE )
2229         cv::boxFilter( src, dst, dst.depth(), cv::Size(param1, param2), cv::Point(-1,-1),
2230             smooth_type == CV_BLUR, cv::BORDER_REPLICATE );
2231     else if( smooth_type == CV_GAUSSIAN )
2232         cv::GaussianBlur( src, dst, cv::Size(param1, param2), param3, param4, cv::BORDER_REPLICATE );
2233     else if( smooth_type == CV_MEDIAN )
2234         cv::medianBlur( src, dst, param1 );
2235     else
2236         cv::bilateralFilter( src, dst, param1, param3, param4, cv::BORDER_REPLICATE );
2237
2238     if( dst.data != dst0.data )
2239         CV_Error( CV_StsUnmatchedFormats, "The destination image does not have the proper type" );
2240 }
2241
2242 /* End of file. */