modules/core/src/stat.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-2011, Willow Garage Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 #include "precomp.hpp"
  44 #include <climits>
  45
  46 namespace cv
  47 {
  48
  49 template<typename T> static inline Scalar rawToScalar(const T& v)
  50 {
  51     Scalar s;
  52     typedef typename DataType<T>::channel_type T1;
  53     int i, n = DataType<T>::channels;
  54     for( i = 0; i < n; i++ )
  55         s.val[i] = ((T1*)&v)[i];
  56     return s;
  57 }
  58
  59 /****************************************************************************************\
  60 *                                        sum                                             *
  61 \****************************************************************************************/
  62
  63 template<typename T, typename ST>
  64 static int sum_(const T* src0, const uchar* mask, ST* dst, int len, int cn )
  65 {
  66     const T* src = src0;
  67     if( !mask )
  68     {
  69         int i=0;
  70         int k = cn % 4;
  71         if( k == 1 )
  72         {
  73             ST s0 = dst[0];
  74
  75             #if CV_ENABLE_UNROLLED
  76             for(; i <= len - 4; i += 4, src += cn*4 )
  77                 s0 += src[0] + src[cn] + src[cn*2] + src[cn*3];
  78             #endif
  79             for( ; i < len; i++, src += cn )
  80                 s0 += src[0];
  81             dst[0] = s0;
  82         }
  83         else if( k == 2 )
  84         {
  85             ST s0 = dst[0], s1 = dst[1];
  86             for( i = 0; i < len; i++, src += cn )
  87             {
  88                 s0 += src[0];
  89                 s1 += src[1];
  90             }
  91             dst[0] = s0;
  92             dst[1] = s1;
  93         }
  94         else if( k == 3 )
  95         {
  96             ST s0 = dst[0], s1 = dst[1], s2 = dst[2];
  97             for( i = 0; i < len; i++, src += cn )
  98             {
  99                 s0 += src[0];
 100                 s1 += src[1];
 101                 s2 += src[2];
 102             }
 103             dst[0] = s0;
 104             dst[1] = s1;
 105             dst[2] = s2;
 106         }
 107
 108         for( ; k < cn; k += 4 )
 109         {
 110             src = src0 + k;
 111             ST s0 = dst[k], s1 = dst[k+1], s2 = dst[k+2], s3 = dst[k+3];
 112             for( i = 0; i < len; i++, src += cn )
 113             {
 114                 s0 += src[0]; s1 += src[1];
 115                 s2 += src[2]; s3 += src[3];
 116             }
 117             dst[k] = s0;
 118             dst[k+1] = s1;
 119             dst[k+2] = s2;
 120             dst[k+3] = s3;
 121         }
 122         return len;
 123     }
 124
 125     int i, nzm = 0;
 126     if( cn == 1 )
 127     {
 128         ST s = dst[0];
 129         for( i = 0; i < len; i++ )
 130             if( mask[i] )
 131             {
 132                 s += src[i];
 133                 nzm++;
 134             }
 135         dst[0] = s;
 136     }
 137     else if( cn == 3 )
 138     {
 139         ST s0 = dst[0], s1 = dst[1], s2 = dst[2];
 140         for( i = 0; i < len; i++, src += 3 )
 141             if( mask[i] )
 142             {
 143                 s0 += src[0];
 144                 s1 += src[1];
 145                 s2 += src[2];
 146                 nzm++;
 147             }
 148         dst[0] = s0;
 149         dst[1] = s1;
 150         dst[2] = s2;
 151     }
 152     else
 153     {
 154         for( i = 0; i < len; i++, src += cn )
 155             if( mask[i] )
 156             {
 157                 int k = 0;
 158                 #if CV_ENABLE_UNROLLED
 159                 for( ; k <= cn - 4; k += 4 )
 160                 {
 161                     ST s0, s1;
 162                     s0 = dst[k] + src[k];
 163                     s1 = dst[k+1] + src[k+1];
 164                     dst[k] = s0; dst[k+1] = s1;
 165                     s0 = dst[k+2] + src[k+2];
 166                     s1 = dst[k+3] + src[k+3];
 167                     dst[k+2] = s0; dst[k+3] = s1;
 168                 }
 169                 #endif
 170                 for( ; k < cn; k++ )
 171                     dst[k] += src[k];
 172                 nzm++;
 173             }
 174     }
 175     return nzm;
 176 }
 177
 178
 179 static int sum8u( const uchar* src, const uchar* mask, int* dst, int len, int cn )
 180 { return sum_(src, mask, dst, len, cn); }
 181
 182 static int sum8s( const schar* src, const uchar* mask, int* dst, int len, int cn )
 183 { return sum_(src, mask, dst, len, cn); }
 184
 185 static int sum16u( const ushort* src, const uchar* mask, int* dst, int len, int cn )
 186 { return sum_(src, mask, dst, len, cn); }
 187
 188 static int sum16s( const short* src, const uchar* mask, int* dst, int len, int cn )
 189 { return sum_(src, mask, dst, len, cn); }
 190
 191 static int sum32s( const int* src, const uchar* mask, double* dst, int len, int cn )
 192 { return sum_(src, mask, dst, len, cn); }
 193
 194 static int sum32f( const float* src, const uchar* mask, double* dst, int len, int cn )
 195 { return sum_(src, mask, dst, len, cn); }
 196
 197 static int sum64f( const double* src, const uchar* mask, double* dst, int len, int cn )
 198 { return sum_(src, mask, dst, len, cn); }
 199
 200 typedef int (*SumFunc)(const uchar*, const uchar* mask, uchar*, int, int);
 201
 202 static SumFunc sumTab[] =
 203 {
 204     (SumFunc)GET_OPTIMIZED(sum8u), (SumFunc)sum8s,
 205     (SumFunc)sum16u, (SumFunc)sum16s,
 206     (SumFunc)sum32s,
 207     (SumFunc)GET_OPTIMIZED(sum32f), (SumFunc)sum64f,
 208     0
 209 };
 210
 211 template<typename T>
 212 static int countNonZero_(const T* src, int len )
 213 {
 214     int i=0, nz = 0;
 215     #if CV_ENABLE_UNROLLED
 216     for(; i <= len - 4; i += 4 )
 217         nz += (src[i] != 0) + (src[i+1] != 0) + (src[i+2] != 0) + (src[i+3] != 0);
 218     #endif
 219     for( ; i < len; i++ )
 220         nz += src[i] != 0;
 221     return nz;
 222 }
 223
 224 static int countNonZero8u( const uchar* src, int len )
 225 {
 226     int i=0, nz = 0;
 227 #if CV_SSE2
 228     if(USE_SSE2)//5x-6x
 229     {
 230         __m128i pattern = _mm_setzero_si128 ();
 231         static uchar tab[256];
 232         static volatile bool initialized = false;
 233         if( !initialized )
 234         {
 235             // we compute inverse popcount table,
 236             // since we pass (img[x] == 0) mask as index in the table.
 237             for( int j = 0; j < 256; j++ )
 238             {
 239                 int val = 0;
 240                 for( int mask = 1; mask < 256; mask += mask )
 241                     val += (j & mask) == 0;
 242                 tab[j] = (uchar)val;
 243             }
 244             initialized = true;
 245         }
 246
 247         for (; i<=len-16; i+=16)
 248         {
 249             __m128i r0 = _mm_loadu_si128((const __m128i*)(src+i));
 250             int val = _mm_movemask_epi8(_mm_cmpeq_epi8(r0, pattern));
 251             nz += tab[val & 255] + tab[val >> 8];
 252         }
 253     }
 254 #endif
 255     for( ; i < len; i++ )
 256         nz += src[i] != 0;
 257     return nz;
 258 }
 259
 260 static int countNonZero16u( const ushort* src, int len )
 261 { return countNonZero_(src, len); }
 262
 263 static int countNonZero32s( const int* src, int len )
 264 { return countNonZero_(src, len); }
 265
 266 static int countNonZero32f( const float* src, int len )
 267 { return countNonZero_(src, len); }
 268
 269 static int countNonZero64f( const double* src, int len )
 270 { return countNonZero_(src, len); }
 271
 272 typedef int (*CountNonZeroFunc)(const uchar*, int);
 273
 274 static CountNonZeroFunc countNonZeroTab[] =
 275 {
 276     (CountNonZeroFunc)GET_OPTIMIZED(countNonZero8u), (CountNonZeroFunc)GET_OPTIMIZED(countNonZero8u),
 277     (CountNonZeroFunc)GET_OPTIMIZED(countNonZero16u), (CountNonZeroFunc)GET_OPTIMIZED(countNonZero16u),
 278     (CountNonZeroFunc)GET_OPTIMIZED(countNonZero32s), (CountNonZeroFunc)GET_OPTIMIZED(countNonZero32f),
 279     (CountNonZeroFunc)GET_OPTIMIZED(countNonZero64f), 0
 280 };
 281
 282
 283 template<typename T, typename ST, typename SQT>
 284 static int sumsqr_(const T* src0, const uchar* mask, ST* sum, SQT* sqsum, int len, int cn )
 285 {
 286     const T* src = src0;
 287
 288     if( !mask )
 289     {
 290         int i;
 291         int k = cn % 4;
 292
 293         if( k == 1 )
 294         {
 295             ST s0 = sum[0];
 296             SQT sq0 = sqsum[0];
 297             for( i = 0; i < len; i++, src += cn )
 298             {
 299                 T v = src[0];
 300                 s0 += v; sq0 += (SQT)v*v;
 301             }
 302             sum[0] = s0;
 303             sqsum[0] = sq0;
 304         }
 305         else if( k == 2 )
 306         {
 307             ST s0 = sum[0], s1 = sum[1];
 308             SQT sq0 = sqsum[0], sq1 = sqsum[1];
 309             for( i = 0; i < len; i++, src += cn )
 310             {
 311                 T v0 = src[0], v1 = src[1];
 312                 s0 += v0; sq0 += (SQT)v0*v0;
 313                 s1 += v1; sq1 += (SQT)v1*v1;
 314             }
 315             sum[0] = s0; sum[1] = s1;
 316             sqsum[0] = sq0; sqsum[1] = sq1;
 317         }
 318         else if( k == 3 )
 319         {
 320             ST s0 = sum[0], s1 = sum[1], s2 = sum[2];
 321             SQT sq0 = sqsum[0], sq1 = sqsum[1], sq2 = sqsum[2];
 322             for( i = 0; i < len; i++, src += cn )
 323             {
 324                 T v0 = src[0], v1 = src[1], v2 = src[2];
 325                 s0 += v0; sq0 += (SQT)v0*v0;
 326                 s1 += v1; sq1 += (SQT)v1*v1;
 327                 s2 += v2; sq2 += (SQT)v2*v2;
 328             }
 329             sum[0] = s0; sum[1] = s1; sum[2] = s2;
 330             sqsum[0] = sq0; sqsum[1] = sq1; sqsum[2] = sq2;
 331         }
 332
 333         for( ; k < cn; k += 4 )
 334         {
 335             src = src0 + k;
 336             ST s0 = sum[k], s1 = sum[k+1], s2 = sum[k+2], s3 = sum[k+3];
 337             SQT sq0 = sqsum[k], sq1 = sqsum[k+1], sq2 = sqsum[k+2], sq3 = sqsum[k+3];
 338             for( i = 0; i < len; i++, src += cn )
 339             {
 340                 T v0, v1;
 341                 v0 = src[0], v1 = src[1];
 342                 s0 += v0; sq0 += (SQT)v0*v0;
 343                 s1 += v1; sq1 += (SQT)v1*v1;
 344                 v0 = src[2], v1 = src[3];
 345                 s2 += v0; sq2 += (SQT)v0*v0;
 346                 s3 += v1; sq3 += (SQT)v1*v1;
 347             }
 348             sum[k] = s0; sum[k+1] = s1;
 349             sum[k+2] = s2; sum[k+3] = s3;
 350             sqsum[k] = sq0; sqsum[k+1] = sq1;
 351             sqsum[k+2] = sq2; sqsum[k+3] = sq3;
 352         }
 353         return len;
 354     }
 355
 356     int i, nzm = 0;
 357
 358     if( cn == 1 )
 359     {
 360         ST s0 = sum[0];
 361         SQT sq0 = sqsum[0];
 362         for( i = 0; i < len; i++ )
 363             if( mask[i] )
 364             {
 365                 T v = src[i];
 366                 s0 += v; sq0 += (SQT)v*v;
 367                 nzm++;
 368             }
 369         sum[0] = s0;
 370         sqsum[0] = sq0;
 371     }
 372     else if( cn == 3 )
 373     {
 374         ST s0 = sum[0], s1 = sum[1], s2 = sum[2];
 375         SQT sq0 = sqsum[0], sq1 = sqsum[1], sq2 = sqsum[2];
 376         for( i = 0; i < len; i++, src += 3 )
 377             if( mask[i] )
 378             {
 379                 T v0 = src[0], v1 = src[1], v2 = src[2];
 380                 s0 += v0; sq0 += (SQT)v0*v0;
 381                 s1 += v1; sq1 += (SQT)v1*v1;
 382                 s2 += v2; sq2 += (SQT)v2*v2;
 383                 nzm++;
 384             }
 385         sum[0] = s0; sum[1] = s1; sum[2] = s2;
 386         sqsum[0] = sq0; sqsum[1] = sq1; sqsum[2] = sq2;
 387     }
 388     else
 389     {
 390         for( i = 0; i < len; i++, src += cn )
 391             if( mask[i] )
 392             {
 393                 for( int k = 0; k < cn; k++ )
 394                 {
 395                     T v = src[k];
 396                     ST s = sum[k] + v;
 397                     SQT sq = sqsum[k] + (SQT)v*v;
 398                     sum[k] = s; sqsum[k] = sq;
 399                 }
 400                 nzm++;
 401             }
 402     }
 403     return nzm;
 404 }
 405
 406
 407 static int sqsum8u( const uchar* src, const uchar* mask, int* sum, int* sqsum, int len, int cn )
 408 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 409
 410 static int sqsum8s( const schar* src, const uchar* mask, int* sum, int* sqsum, int len, int cn )
 411 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 412
 413 static int sqsum16u( const ushort* src, const uchar* mask, int* sum, double* sqsum, int len, int cn )
 414 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 415
 416 static int sqsum16s( const short* src, const uchar* mask, int* sum, double* sqsum, int len, int cn )
 417 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 418
 419 static int sqsum32s( const int* src, const uchar* mask, double* sum, double* sqsum, int len, int cn )
 420 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 421
 422 static int sqsum32f( const float* src, const uchar* mask, double* sum, double* sqsum, int len, int cn )
 423 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 424
 425 static int sqsum64f( const double* src, const uchar* mask, double* sum, double* sqsum, int len, int cn )
 426 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 427
 428 typedef int (*SumSqrFunc)(const uchar*, const uchar* mask, uchar*, uchar*, int, int);
 429
 430 static SumSqrFunc sumSqrTab[] =
 431 {
 432     (SumSqrFunc)GET_OPTIMIZED(sqsum8u), (SumSqrFunc)sqsum8s, (SumSqrFunc)sqsum16u, (SumSqrFunc)sqsum16s,
 433     (SumSqrFunc)sqsum32s, (SumSqrFunc)GET_OPTIMIZED(sqsum32f), (SumSqrFunc)sqsum64f, 0
 434 };
 435
 436 }
 437
 438 cv::Scalar cv::sum( InputArray _src )
 439 {
 440     Mat src = _src.getMat();
 441     int k, cn = src.channels(), depth = src.depth();
 442
 443 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
 444         size_t total_size = src.total();
 445         int rows = src.size[0], cols = (int)(total_size/rows);
 446         if( src.dims == 2 || (src.isContinuous() && cols > 0 && (size_t)rows*cols == total_size) )
 447         {
 448                 IppiSize sz = { cols, rows };
 449                 int type = src.type();
 450                 typedef IppStatus (CV_STDCALL* ippiSumFunc)(const void*, int, IppiSize, double *, int);
 451                 ippiSumFunc ippFunc =
 452                         type == CV_8UC1 ? (ippiSumFunc)ippiSum_8u_C1R :
 453                         type == CV_8UC3 ? (ippiSumFunc)ippiSum_8u_C3R :
 454                         type == CV_8UC4 ? (ippiSumFunc)ippiSum_8u_C4R :
 455                         type == CV_16UC1 ? (ippiSumFunc)ippiSum_16u_C1R :
 456                         type == CV_16UC3 ? (ippiSumFunc)ippiSum_16u_C3R :
 457                         type == CV_16UC4 ? (ippiSumFunc)ippiSum_16u_C4R :
 458                         type == CV_16SC1 ? (ippiSumFunc)ippiSum_16s_C1R :
 459                         type == CV_16SC3 ? (ippiSumFunc)ippiSum_16s_C3R :
 460                         type == CV_16SC4 ? (ippiSumFunc)ippiSum_16s_C4R :
 461                         type == CV_32FC1 ? (ippiSumFunc)ippiSum_32f_C1R :
 462                         type == CV_32FC3 ? (ippiSumFunc)ippiSum_32f_C3R :
 463                         type == CV_32FC4 ? (ippiSumFunc)ippiSum_32f_C4R :
 464                         0;
 465                 if( ippFunc )
 466                 {
 467                         Ipp64f res[4];
 468                         if( ippFunc(src.data, src.step[0], sz, res, ippAlgHintAccurate) >= 0 )
 469                         {
 470                                 Scalar sc;
 471                                 for( int i = 0; i < cn; i++ )
 472                                 {
 473                                         sc[i] = res[i];
 474                                 }
 475                                 return sc;
 476                         }
 477                 }
 478         }
 479 #endif
 480
 481     SumFunc func = sumTab[depth];
 482
 483     CV_Assert( cn <= 4 && func != 0 );
 484
 485     const Mat* arrays[] = {&src, 0};
 486     uchar* ptrs[1];
 487     NAryMatIterator it(arrays, ptrs);
 488     Scalar s;
 489     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
 490     int j, count = 0;
 491     AutoBuffer<int> _buf;
 492     int* buf = (int*)&s[0];
 493     size_t esz = 0;
 494     bool blockSum = depth < CV_32S;
 495
 496     if( blockSum )
 497     {
 498         intSumBlockSize = depth <= CV_8S ? (1 << 23) : (1 << 15);
 499         blockSize = std::min(blockSize, intSumBlockSize);
 500         _buf.allocate(cn);
 501         buf = _buf;
 502
 503         for( k = 0; k < cn; k++ )
 504             buf[k] = 0;
 505         esz = src.elemSize();
 506     }
 507
 508     for( size_t i = 0; i < it.nplanes; i++, ++it )
 509     {
 510         for( j = 0; j < total; j += blockSize )
 511         {
 512             int bsz = std::min(total - j, blockSize);
 513             func( ptrs[0], 0, (uchar*)buf, bsz, cn );
 514             count += bsz;
 515             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
 516             {
 517                 for( k = 0; k < cn; k++ )
 518                 {
 519                     s[k] += buf[k];
 520                     buf[k] = 0;
 521                 }
 522                 count = 0;
 523             }
 524             ptrs[0] += bsz*esz;
 525         }
 526     }
 527     return s;
 528 }
 529
 530 int cv::countNonZero( InputArray _src )
 531 {
 532     Mat src = _src.getMat();
 533     CountNonZeroFunc func = countNonZeroTab[src.depth()];
 534
 535     CV_Assert( src.channels() == 1 && func != 0 );
 536
 537     const Mat* arrays[] = {&src, 0};
 538     uchar* ptrs[1];
 539     NAryMatIterator it(arrays, ptrs);
 540     int total = (int)it.size, nz = 0;
 541
 542     for( size_t i = 0; i < it.nplanes; i++, ++it )
 543         nz += func( ptrs[0], total );
 544
 545     return nz;
 546 }
 547
 548 cv::Scalar cv::mean( InputArray _src, InputArray _mask )
 549 {
 550     Mat src = _src.getMat(), mask = _mask.getMat();
 551     CV_Assert( mask.empty() || mask.type() == CV_8U );
 552
 553     int k, cn = src.channels(), depth = src.depth();
 554
 555 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
 556         size_t total_size = src.total();
 557         int rows = src.size[0], cols = (int)(total_size/rows);
 558         if( src.dims == 2 || (src.isContinuous() && mask.isContinuous() && cols > 0 && (size_t)rows*cols == total_size) )
 559         {
 560                 IppiSize sz = { cols, rows };
 561                 int type = src.type();
 562                 if( !mask.empty() )
 563                 {
 564                         typedef IppStatus (CV_STDCALL* ippiMaskMeanFuncC1)(const void *, int, void *, int, IppiSize, Ipp64f *);
 565                         ippiMaskMeanFuncC1 ippFuncC1 =
 566                         type == CV_8UC1 ? (ippiMaskMeanFuncC1)ippiMean_8u_C1MR :
 567                         type == CV_16UC1 ? (ippiMaskMeanFuncC1)ippiMean_16u_C1MR :
 568                         type == CV_32FC1 ? (ippiMaskMeanFuncC1)ippiMean_32f_C1MR :
 569                         0;
 570                         if( ippFuncC1 )
 571                         {
 572                                 Ipp64f res;
 573                                 if( ippFuncC1(src.data, src.step[0], mask.data, mask.step[0], sz, &res) >= 0 )
 574                                 {
 575                                         return Scalar(res);
 576                                 }
 577                         }
 578                         typedef IppStatus (CV_STDCALL* ippiMaskMeanFuncC3)(const void *, int, void *, int, IppiSize, int, Ipp64f *);
 579                         ippiMaskMeanFuncC3 ippFuncC3 =
 580                         type == CV_8UC3 ? (ippiMaskMeanFuncC3)ippiMean_8u_C3CMR :
 581                         type == CV_16UC3 ? (ippiMaskMeanFuncC3)ippiMean_16u_C3CMR :
 582                         type == CV_32FC3 ? (ippiMaskMeanFuncC3)ippiMean_32f_C3CMR :
 583                         0;
 584                         if( ippFuncC3 )
 585                         {
 586                                 Ipp64f res1, res2, res3;
 587                                 if( ippFuncC3(src.data, src.step[0], mask.data, mask.step[0], sz, 1, &res1) >= 0 &&
 588                                         ippFuncC3(src.data, src.step[0], mask.data, mask.step[0], sz, 2, &res2) >= 0 &&
 589                                         ippFuncC3(src.data, src.step[0], mask.data, mask.step[0], sz, 3, &res3) >= 0 )
 590                                 {
 591                                         return Scalar(res1, res2, res3);
 592                                 }
 593                         }
 594                 }
 595                 else
 596                 {
 597                         typedef IppStatus (CV_STDCALL* ippiMeanFunc)(const void*, int, IppiSize, double *, int);
 598                         ippiMeanFunc ippFunc =
 599                                 type == CV_8UC1 ? (ippiMeanFunc)ippiMean_8u_C1R :
 600                                 type == CV_8UC3 ? (ippiMeanFunc)ippiMean_8u_C3R :
 601                                 type == CV_8UC4 ? (ippiMeanFunc)ippiMean_8u_C4R :
 602                                 type == CV_16UC1 ? (ippiMeanFunc)ippiMean_16u_C1R :
 603                                 type == CV_16UC3 ? (ippiMeanFunc)ippiMean_16u_C3R :
 604                                 type == CV_16UC4 ? (ippiMeanFunc)ippiMean_16u_C4R :
 605                                 type == CV_16SC1 ? (ippiMeanFunc)ippiMean_16s_C1R :
 606                                 type == CV_16SC3 ? (ippiMeanFunc)ippiMean_16s_C3R :
 607                                 type == CV_16SC4 ? (ippiMeanFunc)ippiMean_16s_C4R :
 608                                 type == CV_32FC1 ? (ippiMeanFunc)ippiMean_32f_C1R :
 609                                 type == CV_32FC3 ? (ippiMeanFunc)ippiMean_32f_C3R :
 610                                 type == CV_32FC4 ? (ippiMeanFunc)ippiMean_32f_C4R :
 611                                 0;
 612                         if( ippFunc )
 613                         {
 614                                 Ipp64f res[4];
 615                                 if( ippFunc(src.data, src.step[0], sz, res, ippAlgHintAccurate) >= 0 )
 616                                 {
 617                                         Scalar sc;
 618                                         for( int i = 0; i < cn; i++ )
 619                                         {
 620                                                 sc[i] = res[i];
 621                                         }
 622                                         return sc;
 623                                 }
 624                         }
 625                 }
 626         }
 627 #endif
 628
 629     SumFunc func = sumTab[depth];
 630
 631     CV_Assert( cn <= 4 && func != 0 );
 632
 633     const Mat* arrays[] = {&src, &mask, 0};
 634     uchar* ptrs[2];
 635     NAryMatIterator it(arrays, ptrs);
 636     Scalar s;
 637     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
 638     int j, count = 0;
 639     AutoBuffer<int> _buf;
 640     int* buf = (int*)&s[0];
 641     bool blockSum = depth <= CV_16S;
 642     size_t esz = 0, nz0 = 0;
 643
 644     if( blockSum )
 645     {
 646         intSumBlockSize = depth <= CV_8S ? (1 << 23) : (1 << 15);
 647         blockSize = std::min(blockSize, intSumBlockSize);
 648         _buf.allocate(cn);
 649         buf = _buf;
 650
 651         for( k = 0; k < cn; k++ )
 652             buf[k] = 0;
 653         esz = src.elemSize();
 654     }
 655
 656     for( size_t i = 0; i < it.nplanes; i++, ++it )
 657     {
 658         for( j = 0; j < total; j += blockSize )
 659         {
 660             int bsz = std::min(total - j, blockSize);
 661             int nz = func( ptrs[0], ptrs[1], (uchar*)buf, bsz, cn );
 662             count += nz;
 663             nz0 += nz;
 664             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
 665             {
 666                 for( k = 0; k < cn; k++ )
 667                 {
 668                     s[k] += buf[k];
 669                     buf[k] = 0;
 670                 }
 671                 count = 0;
 672             }
 673             ptrs[0] += bsz*esz;
 674             if( ptrs[1] )
 675                 ptrs[1] += bsz;
 676         }
 677     }
 678     return s*(nz0 ? 1./nz0 : 0);
 679 }
 680
 681
 682 void cv::meanStdDev( InputArray _src, OutputArray _mean, OutputArray _sdv, InputArray _mask )
 683 {
 684     Mat src = _src.getMat(), mask = _mask.getMat();
 685     CV_Assert( mask.empty() || mask.type() == CV_8U );
 686
 687     int k, cn = src.channels(), depth = src.depth();
 688     SumSqrFunc func = sumSqrTab[depth];
 689
 690     CV_Assert( func != 0 );
 691
 692     const Mat* arrays[] = {&src, &mask, 0};
 693     uchar* ptrs[2];
 694     NAryMatIterator it(arrays, ptrs);
 695     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
 696     int j, count = 0, nz0 = 0;
 697     AutoBuffer<double> _buf(cn*4);
 698     double *s = (double*)_buf, *sq = s + cn;
 699     int *sbuf = (int*)s, *sqbuf = (int*)sq;
 700     bool blockSum = depth <= CV_16S, blockSqSum = depth <= CV_8S;
 701     size_t esz = 0;
 702
 703     for( k = 0; k < cn; k++ )
 704         s[k] = sq[k] = 0;
 705
 706     if( blockSum )
 707     {
 708         intSumBlockSize = 1 << 15;
 709         blockSize = std::min(blockSize, intSumBlockSize);
 710         sbuf = (int*)(sq + cn);
 711         if( blockSqSum )
 712             sqbuf = sbuf + cn;
 713         for( k = 0; k < cn; k++ )
 714             sbuf[k] = sqbuf[k] = 0;
 715         esz = src.elemSize();
 716     }
 717
 718     for( size_t i = 0; i < it.nplanes; i++, ++it )
 719     {
 720         for( j = 0; j < total; j += blockSize )
 721         {
 722             int bsz = std::min(total - j, blockSize);
 723             int nz = func( ptrs[0], ptrs[1], (uchar*)sbuf, (uchar*)sqbuf, bsz, cn );
 724             count += nz;
 725             nz0 += nz;
 726             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
 727             {
 728                 for( k = 0; k < cn; k++ )
 729                 {
 730                     s[k] += sbuf[k];
 731                     sbuf[k] = 0;
 732                 }
 733                 if( blockSqSum )
 734                 {
 735                     for( k = 0; k < cn; k++ )
 736                     {
 737                         sq[k] += sqbuf[k];
 738                         sqbuf[k] = 0;
 739                     }
 740                 }
 741                 count = 0;
 742             }
 743             ptrs[0] += bsz*esz;
 744             if( ptrs[1] )
 745                 ptrs[1] += bsz;
 746         }
 747     }
 748
 749     double scale = nz0 ? 1./nz0 : 0.;
 750     for( k = 0; k < cn; k++ )
 751     {
 752         s[k] *= scale;
 753         sq[k] = std::sqrt(std::max(sq[k]*scale - s[k]*s[k], 0.));
 754     }
 755
 756     for( j = 0; j < 2; j++ )
 757     {
 758         const double* sptr = j == 0 ? s : sq;
 759         _OutputArray _dst = j == 0 ? _mean : _sdv;
 760         if( !_dst.needed() )
 761             continue;
 762
 763         if( !_dst.fixedSize() )
 764             _dst.create(cn, 1, CV_64F, -1, true);
 765         Mat dst = _dst.getMat();
 766         int dcn = (int)dst.total();
 767         CV_Assert( dst.type() == CV_64F && dst.isContinuous() &&
 768                    (dst.cols == 1 || dst.rows == 1) && dcn >= cn );
 769         double* dptr = dst.ptr<double>();
 770         for( k = 0; k < cn; k++ )
 771             dptr[k] = sptr[k];
 772         for( ; k < dcn; k++ )
 773             dptr[k] = 0;
 774     }
 775 }
 776
 777 /****************************************************************************************\
 778 *                                       minMaxLoc                                        *
 779 \****************************************************************************************/
 780
 781 namespace cv
 782 {
 783
 784 template<typename T, typename WT> static void
 785 minMaxIdx_( const T* src, const uchar* mask, WT* _minVal, WT* _maxVal,
 786             size_t* _minIdx, size_t* _maxIdx, int len, size_t startIdx )
 787 {
 788     WT minVal = *_minVal, maxVal = *_maxVal;
 789     size_t minIdx = *_minIdx, maxIdx = *_maxIdx;
 790
 791     if( !mask )
 792     {
 793         for( int i = 0; i < len; i++ )
 794         {
 795             T val = src[i];
 796             if( val < minVal )
 797             {
 798                 minVal = val;
 799                 minIdx = startIdx + i;
 800             }
 801             if( val > maxVal )
 802             {
 803                 maxVal = val;
 804                 maxIdx = startIdx + i;
 805             }
 806         }
 807     }
 808     else
 809     {
 810         for( int i = 0; i < len; i++ )
 811         {
 812             T val = src[i];
 813             if( mask[i] && val < minVal )
 814             {
 815                 minVal = val;
 816                 minIdx = startIdx + i;
 817             }
 818             if( mask[i] && val > maxVal )
 819             {
 820                 maxVal = val;
 821                 maxIdx = startIdx + i;
 822             }
 823         }
 824     }
 825
 826     *_minIdx = minIdx;
 827     *_maxIdx = maxIdx;
 828     *_minVal = minVal;
 829     *_maxVal = maxVal;
 830 }
 831
 832 static void minMaxIdx_8u(const uchar* src, const uchar* mask, int* minval, int* maxval,
 833                          size_t* minidx, size_t* maxidx, int len, size_t startidx )
 834 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 835
 836 static void minMaxIdx_8s(const schar* src, const uchar* mask, int* minval, int* maxval,
 837                          size_t* minidx, size_t* maxidx, int len, size_t startidx )
 838 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 839
 840 static void minMaxIdx_16u(const ushort* src, const uchar* mask, int* minval, int* maxval,
 841                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 842 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 843
 844 static void minMaxIdx_16s(const short* src, const uchar* mask, int* minval, int* maxval,
 845                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 846 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 847
 848 static void minMaxIdx_32s(const int* src, const uchar* mask, int* minval, int* maxval,
 849                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 850 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 851
 852 static void minMaxIdx_32f(const float* src, const uchar* mask, float* minval, float* maxval,
 853                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 854 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 855
 856 static void minMaxIdx_64f(const double* src, const uchar* mask, double* minval, double* maxval,
 857                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 858 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 859
 860 typedef void (*MinMaxIdxFunc)(const uchar*, const uchar*, int*, int*, size_t*, size_t*, int, size_t);
 861
 862 static MinMaxIdxFunc minmaxTab[] =
 863 {
 864     (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_8u), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_8s),
 865     (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_16u), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_16s),
 866     (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_32s),
 867     (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_32f), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_64f),
 868     0
 869 };
 870
 871 static void ofs2idx(const Mat& a, size_t ofs, int* idx)
 872 {
 873     int i, d = a.dims;
 874     if( ofs > 0 )
 875     {
 876         ofs--;
 877         for( i = d-1; i >= 0; i-- )
 878         {
 879             int sz = a.size[i];
 880             idx[i] = (int)(ofs % sz);
 881             ofs /= sz;
 882         }
 883     }
 884     else
 885     {
 886         for( i = d-1; i >= 0; i-- )
 887             idx[i] = -1;
 888     }
 889 }
 890
 891 }
 892
 893 void cv::minMaxIdx(InputArray _src, double* minVal,
 894                    double* maxVal, int* minIdx, int* maxIdx,
 895                    InputArray _mask)
 896 {
 897     Mat src = _src.getMat(), mask = _mask.getMat();
 898     int depth = src.depth(), cn = src.channels();
 899
 900     CV_Assert( (cn == 1 && (mask.empty() || mask.type() == CV_8U)) ||
 901                (cn >= 1 && mask.empty() && !minIdx && !maxIdx) );
 902     MinMaxIdxFunc func = minmaxTab[depth];
 903     CV_Assert( func != 0 );
 904
 905     const Mat* arrays[] = {&src, &mask, 0};
 906     uchar* ptrs[2];
 907     NAryMatIterator it(arrays, ptrs);
 908
 909     size_t minidx = 0, maxidx = 0;
 910     int iminval = INT_MAX, imaxval = INT_MIN;
 911     float fminval = FLT_MAX, fmaxval = -FLT_MAX;
 912     double dminval = DBL_MAX, dmaxval = -DBL_MAX;
 913     size_t startidx = 1;
 914     int *minval = &iminval, *maxval = &imaxval;
 915     int planeSize = (int)it.size*cn;
 916
 917     if( depth == CV_32F )
 918         minval = (int*)&fminval, maxval = (int*)&fmaxval;
 919     else if( depth == CV_64F )
 920         minval = (int*)&dminval, maxval = (int*)&dmaxval;
 921
 922     for( size_t i = 0; i < it.nplanes; i++, ++it, startidx += planeSize )
 923         func( ptrs[0], ptrs[1], minval, maxval, &minidx, &maxidx, planeSize, startidx );
 924
 925     if( minidx == 0 )
 926         dminval = dmaxval = 0;
 927     else if( depth == CV_32F )
 928         dminval = fminval, dmaxval = fmaxval;
 929     else if( depth <= CV_32S )
 930         dminval = iminval, dmaxval = imaxval;
 931
 932     if( minVal )
 933         *minVal = dminval;
 934     if( maxVal )
 935         *maxVal = dmaxval;
 936
 937     if( minIdx )
 938         ofs2idx(src, minidx, minIdx);
 939     if( maxIdx )
 940         ofs2idx(src, maxidx, maxIdx);
 941 }
 942
 943 void cv::minMaxLoc( InputArray _img, double* minVal, double* maxVal,
 944                     Point* minLoc, Point* maxLoc, InputArray mask )
 945 {
 946     Mat img = _img.getMat();
 947     CV_Assert(img.dims <= 2);
 948
 949     minMaxIdx(_img, minVal, maxVal, (int*)minLoc, (int*)maxLoc, mask);
 950     if( minLoc )
 951         std::swap(minLoc->x, minLoc->y);
 952     if( maxLoc )
 953         std::swap(maxLoc->x, maxLoc->y);
 954 }
 955
 956 /****************************************************************************************\
 957 *                                         norm                                           *
 958 \****************************************************************************************/
 959
 960 namespace cv
 961 {
 962
 963 float normL2Sqr_(const float* a, const float* b, int n)
 964 {
 965     int j = 0; float d = 0.f;
 966 #if CV_SSE
 967     if( USE_SSE2 )
 968     {
 969         float CV_DECL_ALIGNED(16) buf[4];
 970         __m128 d0 = _mm_setzero_ps(), d1 = _mm_setzero_ps();
 971
 972         for( ; j <= n - 8; j += 8 )
 973         {
 974             __m128 t0 = _mm_sub_ps(_mm_loadu_ps(a + j), _mm_loadu_ps(b + j));
 975             __m128 t1 = _mm_sub_ps(_mm_loadu_ps(a + j + 4), _mm_loadu_ps(b + j + 4));
 976             d0 = _mm_add_ps(d0, _mm_mul_ps(t0, t0));
 977             d1 = _mm_add_ps(d1, _mm_mul_ps(t1, t1));
 978         }
 979         _mm_store_ps(buf, _mm_add_ps(d0, d1));
 980         d = buf[0] + buf[1] + buf[2] + buf[3];
 981     }
 982     else
 983 #endif
 984     {
 985         for( ; j <= n - 4; j += 4 )
 986         {
 987             float t0 = a[j] - b[j], t1 = a[j+1] - b[j+1], t2 = a[j+2] - b[j+2], t3 = a[j+3] - b[j+3];
 988             d += t0*t0 + t1*t1 + t2*t2 + t3*t3;
 989         }
 990     }
 991
 992     for( ; j < n; j++ )
 993     {
 994         float t = a[j] - b[j];
 995         d += t*t;
 996     }
 997     return d;
 998 }
 999
1000
1001 float normL1_(const float* a, const float* b, int n)
1002 {
1003     int j = 0; float d = 0.f;
1004 #if CV_SSE
1005     if( USE_SSE2 )
1006     {
1007         float CV_DECL_ALIGNED(16) buf[4];
1008         static const int CV_DECL_ALIGNED(16) absbuf[4] = {0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff};
1009         __m128 d0 = _mm_setzero_ps(), d1 = _mm_setzero_ps();
1010         __m128 absmask = _mm_load_ps((const float*)absbuf);
1011
1012         for( ; j <= n - 8; j += 8 )
1013         {
1014             __m128 t0 = _mm_sub_ps(_mm_loadu_ps(a + j), _mm_loadu_ps(b + j));
1015             __m128 t1 = _mm_sub_ps(_mm_loadu_ps(a + j + 4), _mm_loadu_ps(b + j + 4));
1016             d0 = _mm_add_ps(d0, _mm_and_ps(t0, absmask));
1017             d1 = _mm_add_ps(d1, _mm_and_ps(t1, absmask));
1018         }
1019         _mm_store_ps(buf, _mm_add_ps(d0, d1));
1020         d = buf[0] + buf[1] + buf[2] + buf[3];
1021     }
1022     else
1023 #endif
1024     {
1025         for( ; j <= n - 4; j += 4 )
1026         {
1027             d += std::abs(a[j] - b[j]) + std::abs(a[j+1] - b[j+1]) +
1028                     std::abs(a[j+2] - b[j+2]) + std::abs(a[j+3] - b[j+3]);
1029         }
1030     }
1031
1032     for( ; j < n; j++ )
1033         d += std::abs(a[j] - b[j]);
1034     return d;
1035 }
1036
1037 int normL1_(const uchar* a, const uchar* b, int n)
1038 {
1039     int j = 0, d = 0;
1040 #if CV_SSE
1041     if( USE_SSE2 )
1042     {
1043         __m128i d0 = _mm_setzero_si128();
1044
1045         for( ; j <= n - 16; j += 16 )
1046         {
1047             __m128i t0 = _mm_loadu_si128((const __m128i*)(a + j));
1048             __m128i t1 = _mm_loadu_si128((const __m128i*)(b + j));
1049
1050             d0 = _mm_add_epi32(d0, _mm_sad_epu8(t0, t1));
1051         }
1052
1053         for( ; j <= n - 4; j += 4 )
1054         {
1055             __m128i t0 = _mm_cvtsi32_si128(*(const int*)(a + j));
1056             __m128i t1 = _mm_cvtsi32_si128(*(const int*)(b + j));
1057
1058             d0 = _mm_add_epi32(d0, _mm_sad_epu8(t0, t1));
1059         }
1060         d = _mm_cvtsi128_si32(_mm_add_epi32(d0, _mm_unpackhi_epi64(d0, d0)));
1061     }
1062     else
1063 #endif
1064     {
1065         for( ; j <= n - 4; j += 4 )
1066         {
1067             d += std::abs(a[j] - b[j]) + std::abs(a[j+1] - b[j+1]) +
1068                     std::abs(a[j+2] - b[j+2]) + std::abs(a[j+3] - b[j+3]);
1069         }
1070     }
1071     for( ; j < n; j++ )
1072         d += std::abs(a[j] - b[j]);
1073     return d;
1074 }
1075
1076 static const uchar popCountTable[] =
1077 {
1078     0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
1079     1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
1080     1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
1081     2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
1082     1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
1083     2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
1084     2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
1085     3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
1086 };
1087
1088 static const uchar popCountTable2[] =
1089 {
1090     0, 1, 1, 1, 1, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3,
1091     1, 2, 2, 2, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 1, 2, 2, 2, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3,
1092     1, 2, 2, 2, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4,
1093     2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4, 2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4,
1094     1, 2, 2, 2, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4,
1095     2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4, 2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4,
1096     1, 2, 2, 2, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4,
1097     2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4, 2, 3, 3, 3, 3, 4, 4, 4, 3, 4, 4, 4, 3, 4, 4, 4
1098 };
1099
1100 static const uchar popCountTable4[] =
1101 {
1102     0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
1103     1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
1104     1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
1105     1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
1106     1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
1107     1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
1108     1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
1109     1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2
1110 };
1111
1112 static int normHamming(const uchar* a, int n)
1113 {
1114     int i = 0, result = 0;
1115 #if CV_NEON
1116     uint32x4_t bits = vmovq_n_u32(0);
1117     for (; i <= n - 16; i += 16) {
1118         uint8x16_t A_vec = vld1q_u8 (a + i);
1119         uint8x16_t bitsSet = vcntq_u8 (A_vec);
1120         uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
1121         uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
1122         bits = vaddq_u32(bits, bitSet4);
1123     }
1124     uint64x2_t bitSet2 = vpaddlq_u32 (bits);
1125     result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
1126     result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
1127 #else
1128     for( ; i <= n - 4; i += 4 )
1129             result += popCountTable[a[i]] + popCountTable[a[i+1]] +
1130             popCountTable[a[i+2]] + popCountTable[a[i+3]];
1131 #endif
1132     for( ; i < n; i++ )
1133         result += popCountTable[a[i]];
1134     return result;
1135 }
1136
1137 int normHamming(const uchar* a, const uchar* b, int n)
1138 {
1139     int i = 0, result = 0;
1140 #if CV_NEON
1141     uint32x4_t bits = vmovq_n_u32(0);
1142     for (; i <= n - 16; i += 16) {
1143         uint8x16_t A_vec = vld1q_u8 (a + i);
1144         uint8x16_t B_vec = vld1q_u8 (b + i);
1145         uint8x16_t AxorB = veorq_u8 (A_vec, B_vec);
1146         uint8x16_t bitsSet = vcntq_u8 (AxorB);
1147         uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
1148         uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
1149         bits = vaddq_u32(bits, bitSet4);
1150     }
1151     uint64x2_t bitSet2 = vpaddlq_u32 (bits);
1152     result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
1153     result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
1154 #else
1155     for( ; i <= n - 4; i += 4 )
1156         result += popCountTable[a[i] ^ b[i]] + popCountTable[a[i+1] ^ b[i+1]] +
1157                 popCountTable[a[i+2] ^ b[i+2]] + popCountTable[a[i+3] ^ b[i+3]];
1158 #endif
1159     for( ; i < n; i++ )
1160         result += popCountTable[a[i] ^ b[i]];
1161     return result;
1162 }
1163
1164 static int normHamming(const uchar* a, int n, int cellSize)
1165 {
1166     if( cellSize == 1 )
1167         return normHamming(a, n);
1168     const uchar* tab = 0;
1169     if( cellSize == 2 )
1170         tab = popCountTable2;
1171     else if( cellSize == 4 )
1172         tab = popCountTable4;
1173     else
1174         CV_Error( CV_StsBadSize, "bad cell size (not 1, 2 or 4) in normHamming" );
1175     int i = 0, result = 0;
1176 #if CV_ENABLE_UNROLLED
1177     for( ; i <= n - 4; i += 4 )
1178         result += tab[a[i]] + tab[a[i+1]] + tab[a[i+2]] + tab[a[i+3]];
1179 #endif
1180     for( ; i < n; i++ )
1181         result += tab[a[i]];
1182     return result;
1183 }
1184
1185 int normHamming(const uchar* a, const uchar* b, int n, int cellSize)
1186 {
1187     if( cellSize == 1 )
1188         return normHamming(a, b, n);
1189     const uchar* tab = 0;
1190     if( cellSize == 2 )
1191         tab = popCountTable2;
1192     else if( cellSize == 4 )
1193         tab = popCountTable4;
1194     else
1195         CV_Error( CV_StsBadSize, "bad cell size (not 1, 2 or 4) in normHamming" );
1196     int i = 0, result = 0;
1197     #if CV_ENABLE_UNROLLED
1198     for( ; i <= n - 4; i += 4 )
1199         result += tab[a[i] ^ b[i]] + tab[a[i+1] ^ b[i+1]] +
1200                 tab[a[i+2] ^ b[i+2]] + tab[a[i+3] ^ b[i+3]];
1201     #endif
1202     for( ; i < n; i++ )
1203         result += tab[a[i] ^ b[i]];
1204     return result;
1205 }
1206
1207
1208 template<typename T, typename ST> int
1209 normInf_(const T* src, const uchar* mask, ST* _result, int len, int cn)
1210 {
1211     ST result = *_result;
1212     if( !mask )
1213     {
1214         result = std::max(result, normInf<T, ST>(src, len*cn));
1215     }
1216     else
1217     {
1218         for( int i = 0; i < len; i++, src += cn )
1219             if( mask[i] )
1220             {
1221                 for( int k = 0; k < cn; k++ )
1222                     result = std::max(result, ST(fast_abs(src[k])));
1223             }
1224     }
1225     *_result = result;
1226     return 0;
1227 }
1228
1229 template<typename T, typename ST> int
1230 normL1_(const T* src, const uchar* mask, ST* _result, int len, int cn)
1231 {
1232     ST result = *_result;
1233     if( !mask )
1234     {
1235         result += normL1<T, ST>(src, len*cn);
1236     }
1237     else
1238     {
1239         for( int i = 0; i < len; i++, src += cn )
1240             if( mask[i] )
1241             {
1242                 for( int k = 0; k < cn; k++ )
1243                     result += fast_abs(src[k]);
1244             }
1245     }
1246     *_result = result;
1247     return 0;
1248 }
1249
1250 template<typename T, typename ST> int
1251 normL2_(const T* src, const uchar* mask, ST* _result, int len, int cn)
1252 {
1253     ST result = *_result;
1254     if( !mask )
1255     {
1256         result += normL2Sqr<T, ST>(src, len*cn);
1257     }
1258     else
1259     {
1260         for( int i = 0; i < len; i++, src += cn )
1261             if( mask[i] )
1262             {
1263                 for( int k = 0; k < cn; k++ )
1264                 {
1265                     T v = src[k];
1266                     result += (ST)v*v;
1267                 }
1268             }
1269     }
1270     *_result = result;
1271     return 0;
1272 }
1273
1274 template<typename T, typename ST> int
1275 normDiffInf_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
1276 {
1277     ST result = *_result;
1278     if( !mask )
1279     {
1280         result = std::max(result, normInf<T, ST>(src1, src2, len*cn));
1281     }
1282     else
1283     {
1284         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
1285             if( mask[i] )
1286             {
1287                 for( int k = 0; k < cn; k++ )
1288                     result = std::max(result, (ST)std::abs(src1[k] - src2[k]));
1289             }
1290     }
1291     *_result = result;
1292     return 0;
1293 }
1294
1295 template<typename T, typename ST> int
1296 normDiffL1_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
1297 {
1298     ST result = *_result;
1299     if( !mask )
1300     {
1301         result += normL1<T, ST>(src1, src2, len*cn);
1302     }
1303     else
1304     {
1305         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
1306             if( mask[i] )
1307             {
1308                 for( int k = 0; k < cn; k++ )
1309                     result += std::abs(src1[k] - src2[k]);
1310             }
1311     }
1312     *_result = result;
1313     return 0;
1314 }
1315
1316 template<typename T, typename ST> int
1317 normDiffL2_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
1318 {
1319     ST result = *_result;
1320     if( !mask )
1321     {
1322         result += normL2Sqr<T, ST>(src1, src2, len*cn);
1323     }
1324     else
1325     {
1326         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
1327             if( mask[i] )
1328             {
1329                 for( int k = 0; k < cn; k++ )
1330                 {
1331                     ST v = src1[k] - src2[k];
1332                     result += v*v;
1333                 }
1334             }
1335     }
1336     *_result = result;
1337     return 0;
1338 }
1339
1340
1341 #define CV_DEF_NORM_FUNC(L, suffix, type, ntype) \
1342     static int norm##L##_##suffix(const type* src, const uchar* mask, ntype* r, int len, int cn) \
1343 { return norm##L##_(src, mask, r, len, cn); } \
1344     static int normDiff##L##_##suffix(const type* src1, const type* src2, \
1345     const uchar* mask, ntype* r, int len, int cn) \
1346 { return normDiff##L##_(src1, src2, mask, r, (int)len, cn); }
1347
1348 #define CV_DEF_NORM_ALL(suffix, type, inftype, l1type, l2type) \
1349     CV_DEF_NORM_FUNC(Inf, suffix, type, inftype) \
1350     CV_DEF_NORM_FUNC(L1, suffix, type, l1type) \
1351     CV_DEF_NORM_FUNC(L2, suffix, type, l2type)
1352
1353 CV_DEF_NORM_ALL(8u, uchar, int, int, int)
1354 CV_DEF_NORM_ALL(8s, schar, int, int, int)
1355 CV_DEF_NORM_ALL(16u, ushort, int, int, double)
1356 CV_DEF_NORM_ALL(16s, short, int, int, double)
1357 CV_DEF_NORM_ALL(32s, int, int, double, double)
1358 CV_DEF_NORM_ALL(32f, float, float, double, double)
1359 CV_DEF_NORM_ALL(64f, double, double, double, double)
1360
1361
1362 typedef int (*NormFunc)(const uchar*, const uchar*, uchar*, int, int);
1363 typedef int (*NormDiffFunc)(const uchar*, const uchar*, const uchar*, uchar*, int, int);
1364
1365 static NormFunc normTab[3][8] =
1366 {
1367     {
1368         (NormFunc)GET_OPTIMIZED(normInf_8u), (NormFunc)GET_OPTIMIZED(normInf_8s), (NormFunc)GET_OPTIMIZED(normInf_16u), (NormFunc)GET_OPTIMIZED(normInf_16s),
1369         (NormFunc)GET_OPTIMIZED(normInf_32s), (NormFunc)GET_OPTIMIZED(normInf_32f), (NormFunc)normInf_64f, 0
1370     },
1371     {
1372         (NormFunc)GET_OPTIMIZED(normL1_8u), (NormFunc)GET_OPTIMIZED(normL1_8s), (NormFunc)GET_OPTIMIZED(normL1_16u), (NormFunc)GET_OPTIMIZED(normL1_16s),
1373         (NormFunc)GET_OPTIMIZED(normL1_32s), (NormFunc)GET_OPTIMIZED(normL1_32f), (NormFunc)normL1_64f, 0
1374     },
1375     {
1376         (NormFunc)GET_OPTIMIZED(normL2_8u), (NormFunc)GET_OPTIMIZED(normL2_8s), (NormFunc)GET_OPTIMIZED(normL2_16u), (NormFunc)GET_OPTIMIZED(normL2_16s),
1377         (NormFunc)GET_OPTIMIZED(normL2_32s), (NormFunc)GET_OPTIMIZED(normL2_32f), (NormFunc)normL2_64f, 0
1378     }
1379 };
1380
1381 static NormDiffFunc normDiffTab[3][8] =
1382 {
1383     {
1384         (NormDiffFunc)GET_OPTIMIZED(normDiffInf_8u), (NormDiffFunc)normDiffInf_8s,
1385         (NormDiffFunc)normDiffInf_16u, (NormDiffFunc)normDiffInf_16s,
1386         (NormDiffFunc)normDiffInf_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffInf_32f),
1387         (NormDiffFunc)normDiffInf_64f, 0
1388     },
1389     {
1390         (NormDiffFunc)GET_OPTIMIZED(normDiffL1_8u), (NormDiffFunc)normDiffL1_8s,
1391         (NormDiffFunc)normDiffL1_16u, (NormDiffFunc)normDiffL1_16s,
1392         (NormDiffFunc)normDiffL1_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffL1_32f),
1393         (NormDiffFunc)normDiffL1_64f, 0
1394     },
1395     {
1396         (NormDiffFunc)GET_OPTIMIZED(normDiffL2_8u), (NormDiffFunc)normDiffL2_8s,
1397         (NormDiffFunc)normDiffL2_16u, (NormDiffFunc)normDiffL2_16s,
1398         (NormDiffFunc)normDiffL2_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffL2_32f),
1399         (NormDiffFunc)normDiffL2_64f, 0
1400     }
1401 };
1402
1403 }
1404
1405 double cv::norm( InputArray _src, int normType, InputArray _mask )
1406 {
1407     Mat src = _src.getMat(), mask = _mask.getMat();
1408     int depth = src.depth(), cn = src.channels();
1409
1410     normType &= 7;
1411     CV_Assert( normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 || normType == NORM_L2SQR ||
1412                ((normType == NORM_HAMMING || normType == NORM_HAMMING2) && src.type() == CV_8U) );
1413
1414     if( src.isContinuous() && mask.empty() )
1415     {
1416         size_t len = src.total()*cn;
1417         if( len == (size_t)(int)len )
1418         {
1419             if( depth == CV_32F )
1420             {
1421                 const float* data = src.ptr<float>();
1422
1423                 if( normType == NORM_L2 )
1424                 {
1425                     double result = 0;
1426                     GET_OPTIMIZED(normL2_32f)(data, 0, &result, (int)len, 1);
1427                     return std::sqrt(result);
1428                 }
1429                 if( normType == NORM_L2SQR )
1430                 {
1431                     double result = 0;
1432                     GET_OPTIMIZED(normL2_32f)(data, 0, &result, (int)len, 1);
1433                     return result;
1434                 }
1435                 if( normType == NORM_L1 )
1436                 {
1437                     double result = 0;
1438                     GET_OPTIMIZED(normL1_32f)(data, 0, &result, (int)len, 1);
1439                     return result;
1440                 }
1441                 if( normType == NORM_INF )
1442                 {
1443                     float result = 0;
1444                     GET_OPTIMIZED(normInf_32f)(data, 0, &result, (int)len, 1);
1445                     return result;
1446                 }
1447             }
1448             if( depth == CV_8U )
1449             {
1450                 const uchar* data = src.ptr<uchar>();
1451
1452                 if( normType == NORM_HAMMING )
1453                     return normHamming(data, (int)len);
1454
1455                 if( normType == NORM_HAMMING2 )
1456                     return normHamming(data, (int)len, 2);
1457             }
1458         }
1459     }
1460
1461     CV_Assert( mask.empty() || mask.type() == CV_8U );
1462
1463     if( normType == NORM_HAMMING || normType == NORM_HAMMING2 )
1464     {
1465         if( !mask.empty() )
1466         {
1467             Mat temp;
1468             bitwise_and(src, mask, temp);
1469             return norm(temp, normType);
1470         }
1471         int cellSize = normType == NORM_HAMMING ? 1 : 2;
1472
1473         const Mat* arrays[] = {&src, 0};
1474         uchar* ptrs[1];
1475         NAryMatIterator it(arrays, ptrs);
1476         int total = (int)it.size;
1477         int result = 0;
1478
1479         for( size_t i = 0; i < it.nplanes; i++, ++it )
1480             result += normHamming(ptrs[0], total, cellSize);
1481
1482         return result;
1483     }
1484
1485     NormFunc func = normTab[normType >> 1][depth];
1486     CV_Assert( func != 0 );
1487
1488     const Mat* arrays[] = {&src, &mask, 0};
1489     uchar* ptrs[2];
1490     union
1491     {
1492         double d;
1493         int i;
1494         float f;
1495     }
1496     result;
1497     result.d = 0;
1498     NAryMatIterator it(arrays, ptrs);
1499     int j, total = (int)it.size, blockSize = total, intSumBlockSize = 0, count = 0;
1500     bool blockSum = (normType == NORM_L1 && depth <= CV_16S) ||
1501             ((normType == NORM_L2 || normType == NORM_L2SQR) && depth <= CV_8S);
1502     int isum = 0;
1503     int *ibuf = &result.i;
1504     size_t esz = 0;
1505
1506     if( blockSum )
1507     {
1508         intSumBlockSize = (normType == NORM_L1 && depth <= CV_8S ? (1 << 23) : (1 << 15))/cn;
1509         blockSize = std::min(blockSize, intSumBlockSize);
1510         ibuf = &isum;
1511         esz = src.elemSize();
1512     }
1513
1514     for( size_t i = 0; i < it.nplanes; i++, ++it )
1515     {
1516         for( j = 0; j < total; j += blockSize )
1517         {
1518             int bsz = std::min(total - j, blockSize);
1519             func( ptrs[0], ptrs[1], (uchar*)ibuf, bsz, cn );
1520             count += bsz;
1521             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
1522             {
1523                 result.d += isum;
1524                 isum = 0;
1525                 count = 0;
1526             }
1527             ptrs[0] += bsz*esz;
1528             if( ptrs[1] )
1529                 ptrs[1] += bsz;
1530         }
1531     }
1532
1533     if( normType == NORM_INF )
1534     {
1535         if( depth == CV_64F )
1536             ;
1537         else if( depth == CV_32F )
1538             result.d = result.f;
1539         else
1540             result.d = result.i;
1541     }
1542     else if( normType == NORM_L2 )
1543         result.d = std::sqrt(result.d);
1544
1545     return result.d;
1546 }
1547
1548
1549 double cv::norm( InputArray _src1, InputArray _src2, int normType, InputArray _mask )
1550 {
1551     if( normType & CV_RELATIVE )
1552         return norm(_src1, _src2, normType & ~CV_RELATIVE, _mask)/(norm(_src2, normType, _mask) + DBL_EPSILON);
1553
1554     Mat src1 = _src1.getMat(), src2 = _src2.getMat(), mask = _mask.getMat();
1555     int depth = src1.depth(), cn = src1.channels();
1556
1557     CV_Assert( src1.size == src2.size && src1.type() == src2.type() );
1558
1559     normType &= 7;
1560     CV_Assert( normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 || normType == NORM_L2SQR ||
1561               ((normType == NORM_HAMMING || normType == NORM_HAMMING2) && src1.type() == CV_8U) );
1562
1563     if( src1.isContinuous() && src2.isContinuous() && mask.empty() )
1564     {
1565         size_t len = src1.total()*src1.channels();
1566         if( len == (size_t)(int)len )
1567         {
1568             if( src1.depth() == CV_32F )
1569             {
1570                 const float* data1 = src1.ptr<float>();
1571                 const float* data2 = src2.ptr<float>();
1572
1573                 if( normType == NORM_L2 )
1574                 {
1575                     double result = 0;
1576                     GET_OPTIMIZED(normDiffL2_32f)(data1, data2, 0, &result, (int)len, 1);
1577                     return std::sqrt(result);
1578                 }
1579                 if( normType == NORM_L2SQR )
1580                 {
1581                     double result = 0;
1582                     GET_OPTIMIZED(normDiffL2_32f)(data1, data2, 0, &result, (int)len, 1);
1583                     return result;
1584                 }
1585                 if( normType == NORM_L1 )
1586                 {
1587                     double result = 0;
1588                     GET_OPTIMIZED(normDiffL1_32f)(data1, data2, 0, &result, (int)len, 1);
1589                     return result;
1590                 }
1591                 if( normType == NORM_INF )
1592                 {
1593                     float result = 0;
1594                     GET_OPTIMIZED(normDiffInf_32f)(data1, data2, 0, &result, (int)len, 1);
1595                     return result;
1596                 }
1597             }
1598         }
1599     }
1600
1601     CV_Assert( mask.empty() || mask.type() == CV_8U );
1602
1603     if( normType == NORM_HAMMING || normType == NORM_HAMMING2 )
1604     {
1605         if( !mask.empty() )
1606         {
1607             Mat temp;
1608             bitwise_xor(src1, src2, temp);
1609             bitwise_and(temp, mask, temp);
1610             return norm(temp, normType);
1611         }
1612         int cellSize = normType == NORM_HAMMING ? 1 : 2;
1613
1614         const Mat* arrays[] = {&src1, &src2, 0};
1615         uchar* ptrs[2];
1616         NAryMatIterator it(arrays, ptrs);
1617         int total = (int)it.size;
1618         int result = 0;
1619
1620         for( size_t i = 0; i < it.nplanes; i++, ++it )
1621             result += normHamming(ptrs[0], ptrs[1], total, cellSize);
1622
1623         return result;
1624     }
1625
1626     NormDiffFunc func = normDiffTab[normType >> 1][depth];
1627     CV_Assert( func != 0 );
1628
1629     const Mat* arrays[] = {&src1, &src2, &mask, 0};
1630     uchar* ptrs[3];
1631     union
1632     {
1633         double d;
1634         float f;
1635         int i;
1636         unsigned u;
1637     }
1638     result;
1639     result.d = 0;
1640     NAryMatIterator it(arrays, ptrs);
1641     int j, total = (int)it.size, blockSize = total, intSumBlockSize = 0, count = 0;
1642     bool blockSum = (normType == NORM_L1 && depth <= CV_16S) ||
1643             ((normType == NORM_L2 || normType == NORM_L2SQR) && depth <= CV_8S);
1644     unsigned isum = 0;
1645     unsigned *ibuf = &result.u;
1646     size_t esz = 0;
1647
1648     if( blockSum )
1649     {
1650         intSumBlockSize = normType == NORM_L1 && depth <= CV_8S ? (1 << 23) : (1 << 15);
1651         blockSize = std::min(blockSize, intSumBlockSize);
1652         ibuf = &isum;
1653         esz = src1.elemSize();
1654     }
1655
1656     for( size_t i = 0; i < it.nplanes; i++, ++it )
1657     {
1658         for( j = 0; j < total; j += blockSize )
1659         {
1660             int bsz = std::min(total - j, blockSize);
1661             func( ptrs[0], ptrs[1], ptrs[2], (uchar*)ibuf, bsz, cn );
1662             count += bsz;
1663             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
1664             {
1665                 result.d += isum;
1666                 isum = 0;
1667                 count = 0;
1668             }
1669             ptrs[0] += bsz*esz;
1670             ptrs[1] += bsz*esz;
1671             if( ptrs[2] )
1672                 ptrs[2] += bsz;
1673         }
1674     }
1675
1676     if( normType == NORM_INF )
1677     {
1678         if( depth == CV_64F )
1679             ;
1680         else if( depth == CV_32F )
1681             result.d = result.f;
1682         else
1683             result.d = result.u;
1684     }
1685     else if( normType == NORM_L2 )
1686         result.d = std::sqrt(result.d);
1687
1688     return result.d;
1689 }
1690
1691
1692 ///////////////////////////////////// batch distance ///////////////////////////////////////
1693
1694 namespace cv
1695 {
1696
1697 template<typename _Tp, typename _Rt>
1698 void batchDistL1_(const _Tp* src1, const _Tp* src2, size_t step2,
1699                   int nvecs, int len, _Rt* dist, const uchar* mask)
1700 {
1701     step2 /= sizeof(src2[0]);
1702     if( !mask )
1703     {
1704         for( int i = 0; i < nvecs; i++ )
1705             dist[i] = normL1<_Tp, _Rt>(src1, src2 + step2*i, len);
1706     }
1707     else
1708     {
1709         _Rt val0 = std::numeric_limits<_Rt>::max();
1710         for( int i = 0; i < nvecs; i++ )
1711             dist[i] = mask[i] ? normL1<_Tp, _Rt>(src1, src2 + step2*i, len) : val0;
1712     }
1713 }
1714
1715 template<typename _Tp, typename _Rt>
1716 void batchDistL2Sqr_(const _Tp* src1, const _Tp* src2, size_t step2,
1717                      int nvecs, int len, _Rt* dist, const uchar* mask)
1718 {
1719     step2 /= sizeof(src2[0]);
1720     if( !mask )
1721     {
1722         for( int i = 0; i < nvecs; i++ )
1723             dist[i] = normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len);
1724     }
1725     else
1726     {
1727         _Rt val0 = std::numeric_limits<_Rt>::max();
1728         for( int i = 0; i < nvecs; i++ )
1729             dist[i] = mask[i] ? normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len) : val0;
1730     }
1731 }
1732
1733 template<typename _Tp, typename _Rt>
1734 void batchDistL2_(const _Tp* src1, const _Tp* src2, size_t step2,
1735                   int nvecs, int len, _Rt* dist, const uchar* mask)
1736 {
1737     step2 /= sizeof(src2[0]);
1738     if( !mask )
1739     {
1740         for( int i = 0; i < nvecs; i++ )
1741             dist[i] = std::sqrt(normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len));
1742     }
1743     else
1744     {
1745         _Rt val0 = std::numeric_limits<_Rt>::max();
1746         for( int i = 0; i < nvecs; i++ )
1747             dist[i] = mask[i] ? std::sqrt(normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len)) : val0;
1748     }
1749 }
1750
1751 static void batchDistHamming(const uchar* src1, const uchar* src2, size_t step2,
1752                              int nvecs, int len, int* dist, const uchar* mask)
1753 {
1754     step2 /= sizeof(src2[0]);
1755     if( !mask )
1756     {
1757         for( int i = 0; i < nvecs; i++ )
1758             dist[i] = normHamming(src1, src2 + step2*i, len);
1759     }
1760     else
1761     {
1762         int val0 = INT_MAX;
1763         for( int i = 0; i < nvecs; i++ )
1764             dist[i] = mask[i] ? normHamming(src1, src2 + step2*i, len) : val0;
1765     }
1766 }
1767
1768 static void batchDistHamming2(const uchar* src1, const uchar* src2, size_t step2,
1769                               int nvecs, int len, int* dist, const uchar* mask)
1770 {
1771     step2 /= sizeof(src2[0]);
1772     if( !mask )
1773     {
1774         for( int i = 0; i < nvecs; i++ )
1775             dist[i] = normHamming(src1, src2 + step2*i, len, 2);
1776     }
1777     else
1778     {
1779         int val0 = INT_MAX;
1780         for( int i = 0; i < nvecs; i++ )
1781             dist[i] = mask[i] ? normHamming(src1, src2 + step2*i, len, 2) : val0;
1782     }
1783 }
1784
1785 static void batchDistL1_8u32s(const uchar* src1, const uchar* src2, size_t step2,
1786                                int nvecs, int len, int* dist, const uchar* mask)
1787 {
1788     batchDistL1_<uchar, int>(src1, src2, step2, nvecs, len, dist, mask);
1789 }
1790
1791 static void batchDistL1_8u32f(const uchar* src1, const uchar* src2, size_t step2,
1792                                int nvecs, int len, float* dist, const uchar* mask)
1793 {
1794     batchDistL1_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
1795 }
1796
1797 static void batchDistL2Sqr_8u32s(const uchar* src1, const uchar* src2, size_t step2,
1798                                   int nvecs, int len, int* dist, const uchar* mask)
1799 {
1800     batchDistL2Sqr_<uchar, int>(src1, src2, step2, nvecs, len, dist, mask);
1801 }
1802
1803 static void batchDistL2Sqr_8u32f(const uchar* src1, const uchar* src2, size_t step2,
1804                                   int nvecs, int len, float* dist, const uchar* mask)
1805 {
1806     batchDistL2Sqr_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
1807 }
1808
1809 static void batchDistL2_8u32f(const uchar* src1, const uchar* src2, size_t step2,
1810                                int nvecs, int len, float* dist, const uchar* mask)
1811 {
1812     batchDistL2_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
1813 }
1814
1815 static void batchDistL1_32f(const float* src1, const float* src2, size_t step2,
1816                              int nvecs, int len, float* dist, const uchar* mask)
1817 {
1818     batchDistL1_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
1819 }
1820
1821 static void batchDistL2Sqr_32f(const float* src1, const float* src2, size_t step2,
1822                                 int nvecs, int len, float* dist, const uchar* mask)
1823 {
1824     batchDistL2Sqr_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
1825 }
1826
1827 static void batchDistL2_32f(const float* src1, const float* src2, size_t step2,
1828                              int nvecs, int len, float* dist, const uchar* mask)
1829 {
1830     batchDistL2_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
1831 }
1832
1833 typedef void (*BatchDistFunc)(const uchar* src1, const uchar* src2, size_t step2,
1834                               int nvecs, int len, uchar* dist, const uchar* mask);
1835
1836
1837 struct BatchDistInvoker : public ParallelLoopBody
1838 {
1839     BatchDistInvoker( const Mat& _src1, const Mat& _src2,
1840                       Mat& _dist, Mat& _nidx, int _K,
1841                       const Mat& _mask, int _update,
1842                       BatchDistFunc _func)
1843     {
1844         src1 = &_src1;
1845         src2 = &_src2;
1846         dist = &_dist;
1847         nidx = &_nidx;
1848         K = _K;
1849         mask = &_mask;
1850         update = _update;
1851         func = _func;
1852     }
1853
1854     void operator()(const Range& range) const
1855     {
1856         AutoBuffer<int> buf(src2->rows);
1857         int* bufptr = buf;
1858
1859         for( int i = range.start; i < range.end; i++ )
1860         {
1861             func(src1->ptr(i), src2->ptr(), src2->step, src2->rows, src2->cols,
1862                  K > 0 ? (uchar*)bufptr : dist->ptr(i), mask->data ? mask->ptr(i) : 0);
1863
1864             if( K > 0 )
1865             {
1866                 int* nidxptr = nidx->ptr<int>(i);
1867                 // since positive float's can be compared just like int's,
1868                 // we handle both CV_32S and CV_32F cases with a single branch
1869                 int* distptr = (int*)dist->ptr(i);
1870
1871                 int j, k;
1872
1873                 for( j = 0; j < src2->rows; j++ )
1874                 {
1875                     int d = bufptr[j];
1876                     if( d < distptr[K-1] )
1877                     {
1878                         for( k = K-2; k >= 0 && distptr[k] > d; k-- )
1879                         {
1880                             nidxptr[k+1] = nidxptr[k];
1881                             distptr[k+1] = distptr[k];
1882                         }
1883                         nidxptr[k+1] = j + update;
1884                         distptr[k+1] = d;
1885                     }
1886                 }
1887             }
1888         }
1889     }
1890
1891     const Mat *src1;
1892     const Mat *src2;
1893     Mat *dist;
1894     Mat *nidx;
1895     const Mat *mask;
1896     int K;
1897     int update;
1898     BatchDistFunc func;
1899 };
1900
1901 }
1902
1903 void cv::batchDistance( InputArray _src1, InputArray _src2,
1904                         OutputArray _dist, int dtype, OutputArray _nidx,
1905                         int normType, int K, InputArray _mask,
1906                         int update, bool crosscheck )
1907 {
1908     Mat src1 = _src1.getMat(), src2 = _src2.getMat(), mask = _mask.getMat();
1909     int type = src1.type();
1910     CV_Assert( type == src2.type() && src1.cols == src2.cols &&
1911                (type == CV_32F || type == CV_8U));
1912     CV_Assert( _nidx.needed() == (K > 0) );
1913
1914     if( dtype == -1 )
1915     {
1916         dtype = normType == NORM_HAMMING || normType == NORM_HAMMING2 ? CV_32S : CV_32F;
1917     }
1918     CV_Assert( (type == CV_8U && dtype == CV_32S) || dtype == CV_32F);
1919
1920     K = std::min(K, src2.rows);
1921
1922     _dist.create(src1.rows, (K > 0 ? K : src2.rows), dtype);
1923     Mat dist = _dist.getMat(), nidx;
1924     if( _nidx.needed() )
1925     {
1926         _nidx.create(dist.size(), CV_32S);
1927         nidx = _nidx.getMat();
1928     }
1929
1930     if( update == 0 && K > 0 )
1931     {
1932         dist = Scalar::all(dtype == CV_32S ? (double)INT_MAX : (double)FLT_MAX);
1933         nidx = Scalar::all(-1);
1934     }
1935
1936     if( crosscheck )
1937     {
1938         CV_Assert( K == 1 && update == 0 && mask.empty() );
1939         Mat tdist, tidx;
1940         batchDistance(src2, src1, tdist, dtype, tidx, normType, K, mask, 0, false);
1941
1942         // if an idx-th element from src1 appeared to be the nearest to i-th element of src2,
1943         // we update the minimum mutual distance between idx-th element of src1 and the whole src2 set.
1944         // As a result, if nidx[idx] = i*, it means that idx-th element of src1 is the nearest
1945         // to i*-th element of src2 and i*-th element of src2 is the closest to idx-th element of src1.
1946         // If nidx[idx] = -1, it means that there is no such ideal couple for it in src2.
1947         // This O(N) procedure is called cross-check and it helps to eliminate some false matches.
1948         if( dtype == CV_32S )
1949         {
1950             for( int i = 0; i < tdist.rows; i++ )
1951             {
1952                 int idx = tidx.at<int>(i);
1953                 int d = tdist.at<int>(i), d0 = dist.at<int>(idx);
1954                 if( d < d0 )
1955                 {
1956                     dist.at<int>(idx) = d;
1957                     nidx.at<int>(idx) = i + update;
1958                 }
1959             }
1960         }
1961         else
1962         {
1963             for( int i = 0; i < tdist.rows; i++ )
1964             {
1965                 int idx = tidx.at<int>(i);
1966                 float d = tdist.at<float>(i), d0 = dist.at<float>(idx);
1967                 if( d < d0 )
1968                 {
1969                     dist.at<float>(idx) = d;
1970                     nidx.at<int>(idx) = i + update;
1971                 }
1972             }
1973         }
1974         return;
1975     }
1976
1977     BatchDistFunc func = 0;
1978     if( type == CV_8U )
1979     {
1980         if( normType == NORM_L1 && dtype == CV_32S )
1981             func = (BatchDistFunc)batchDistL1_8u32s;
1982         else if( normType == NORM_L1 && dtype == CV_32F )
1983             func = (BatchDistFunc)batchDistL1_8u32f;
1984         else if( normType == NORM_L2SQR && dtype == CV_32S )
1985             func = (BatchDistFunc)batchDistL2Sqr_8u32s;
1986         else if( normType == NORM_L2SQR && dtype == CV_32F )
1987             func = (BatchDistFunc)batchDistL2Sqr_8u32f;
1988         else if( normType == NORM_L2 && dtype == CV_32F )
1989             func = (BatchDistFunc)batchDistL2_8u32f;
1990         else if( normType == NORM_HAMMING && dtype == CV_32S )
1991             func = (BatchDistFunc)batchDistHamming;
1992         else if( normType == NORM_HAMMING2 && dtype == CV_32S )
1993             func = (BatchDistFunc)batchDistHamming2;
1994     }
1995     else if( type == CV_32F && dtype == CV_32F )
1996     {
1997         if( normType == NORM_L1 )
1998             func = (BatchDistFunc)batchDistL1_32f;
1999         else if( normType == NORM_L2SQR )
2000             func = (BatchDistFunc)batchDistL2Sqr_32f;
2001         else if( normType == NORM_L2 )
2002             func = (BatchDistFunc)batchDistL2_32f;
2003     }
2004
2005     if( func == 0 )
2006         CV_Error_(CV_StsUnsupportedFormat,
2007                   ("The combination of type=%d, dtype=%d and normType=%d is not supported",
2008                    type, dtype, normType));
2009
2010     parallel_for_(Range(0, src1.rows),
2011                   BatchDistInvoker(src1, src2, dist, nidx, K, mask, update, func));
2012 }
2013
2014
2015 void cv::findNonZero( InputArray _src, OutputArray _idx )
2016 {
2017     Mat src = _src.getMat();
2018     CV_Assert( src.type() == CV_8UC1 );
2019     int n = countNonZero(src);
2020     if( _idx.kind() == _InputArray::MAT && !_idx.getMatRef().isContinuous() )
2021         _idx.release();
2022     _idx.create(n, 1, CV_32SC2);
2023     Mat idx = _idx.getMat();
2024     CV_Assert(idx.isContinuous());
2025     Point* idx_ptr = (Point*)idx.data;
2026
2027     for( int i = 0; i < src.rows; i++ )
2028     {
2029         const uchar* bin_ptr = src.ptr(i);
2030         for( int j = 0; j < src.cols; j++ )
2031             if( bin_ptr[j] )
2032                 *idx_ptr++ = Point(j, i);
2033     }
2034 }
2035
2036
2037 CV_IMPL CvScalar cvSum( const CvArr* srcarr )
2038 {
2039     cv::Scalar sum = cv::sum(cv::cvarrToMat(srcarr, false, true, 1));
2040     if( CV_IS_IMAGE(srcarr) )
2041     {
2042         int coi = cvGetImageCOI((IplImage*)srcarr);
2043         if( coi )
2044         {
2045             CV_Assert( 0 < coi && coi <= 4 );
2046             sum = cv::Scalar(sum[coi-1]);
2047         }
2048     }
2049     return sum;
2050 }
2051
2052 CV_IMPL int cvCountNonZero( const CvArr* imgarr )
2053 {
2054     cv::Mat img = cv::cvarrToMat(imgarr, false, true, 1);
2055     if( img.channels() > 1 )
2056         cv::extractImageCOI(imgarr, img);
2057     return countNonZero(img);
2058 }
2059
2060
2061 CV_IMPL  CvScalar
2062 cvAvg( const void* imgarr, const void* maskarr )
2063 {
2064     cv::Mat img = cv::cvarrToMat(imgarr, false, true, 1);
2065     cv::Scalar mean = !maskarr ? cv::mean(img) : cv::mean(img, cv::cvarrToMat(maskarr));
2066     if( CV_IS_IMAGE(imgarr) )
2067     {
2068         int coi = cvGetImageCOI((IplImage*)imgarr);
2069         if( coi )
2070         {
2071             CV_Assert( 0 < coi && coi <= 4 );
2072             mean = cv::Scalar(mean[coi-1]);
2073         }
2074     }
2075     return mean;
2076 }
2077
2078
2079 CV_IMPL  void
2080 cvAvgSdv( const CvArr* imgarr, CvScalar* _mean, CvScalar* _sdv, const void* maskarr )
2081 {
2082     cv::Scalar mean, sdv;
2083
2084     cv::Mat mask;
2085     if( maskarr )
2086         mask = cv::cvarrToMat(maskarr);
2087
2088     cv::meanStdDev(cv::cvarrToMat(imgarr, false, true, 1), mean, sdv, mask );
2089
2090     if( CV_IS_IMAGE(imgarr) )
2091     {
2092         int coi = cvGetImageCOI((IplImage*)imgarr);
2093         if( coi )
2094         {
2095             CV_Assert( 0 < coi && coi <= 4 );
2096             mean = cv::Scalar(mean[coi-1]);
2097             sdv = cv::Scalar(sdv[coi-1]);
2098         }
2099     }
2100
2101     if( _mean )
2102         *(cv::Scalar*)_mean = mean;
2103     if( _sdv )
2104         *(cv::Scalar*)_sdv = sdv;
2105 }
2106
2107
2108 CV_IMPL void
2109 cvMinMaxLoc( const void* imgarr, double* _minVal, double* _maxVal,
2110              CvPoint* _minLoc, CvPoint* _maxLoc, const void* maskarr )
2111 {
2112     cv::Mat mask, img = cv::cvarrToMat(imgarr, false, true, 1);
2113     if( maskarr )
2114         mask = cv::cvarrToMat(maskarr);
2115     if( img.channels() > 1 )
2116         cv::extractImageCOI(imgarr, img);
2117
2118     cv::minMaxLoc( img, _minVal, _maxVal,
2119                    (cv::Point*)_minLoc, (cv::Point*)_maxLoc, mask );
2120 }
2121
2122
2123 CV_IMPL  double
2124 cvNorm( const void* imgA, const void* imgB, int normType, const void* maskarr )
2125 {
2126     cv::Mat a, mask;
2127     if( !imgA )
2128     {
2129         imgA = imgB;
2130         imgB = 0;
2131     }
2132
2133     a = cv::cvarrToMat(imgA, false, true, 1);
2134     if( maskarr )
2135         mask = cv::cvarrToMat(maskarr);
2136
2137     if( a.channels() > 1 && CV_IS_IMAGE(imgA) && cvGetImageCOI((const IplImage*)imgA) > 0 )
2138         cv::extractImageCOI(imgA, a);
2139
2140     if( !imgB )
2141         return !maskarr ? cv::norm(a, normType) : cv::norm(a, normType, mask);
2142
2143     cv::Mat b = cv::cvarrToMat(imgB, false, true, 1);
2144     if( b.channels() > 1 && CV_IS_IMAGE(imgB) && cvGetImageCOI((const IplImage*)imgB) > 0 )
2145         cv::extractImageCOI(imgB, b);
2146
2147     return !maskarr ? cv::norm(a, b, normType) : cv::norm(a, b, normType, mask);
2148 }