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     SumFunc func = sumTab[depth];
 443
 444     CV_Assert( cn <= 4 && func != 0 );
 445
 446     const Mat* arrays[] = {&src, 0};
 447     uchar* ptrs[1];
 448     NAryMatIterator it(arrays, ptrs);
 449     Scalar s;
 450     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
 451     int j, count = 0;
 452     AutoBuffer<int> _buf;
 453     int* buf = (int*)&s[0];
 454     size_t esz = 0;
 455     bool blockSum = depth < CV_32S;
 456
 457     if( blockSum )
 458     {
 459         intSumBlockSize = depth <= CV_8S ? (1 << 23) : (1 << 15);
 460         blockSize = std::min(blockSize, intSumBlockSize);
 461         _buf.allocate(cn);
 462         buf = _buf;
 463
 464         for( k = 0; k < cn; k++ )
 465             buf[k] = 0;
 466         esz = src.elemSize();
 467     }
 468
 469     for( size_t i = 0; i < it.nplanes; i++, ++it )
 470     {
 471         for( j = 0; j < total; j += blockSize )
 472         {
 473             int bsz = std::min(total - j, blockSize);
 474             func( ptrs[0], 0, (uchar*)buf, bsz, cn );
 475             count += bsz;
 476             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
 477             {
 478                 for( k = 0; k < cn; k++ )
 479                 {
 480                     s[k] += buf[k];
 481                     buf[k] = 0;
 482                 }
 483                 count = 0;
 484             }
 485             ptrs[0] += bsz*esz;
 486         }
 487     }
 488     return s;
 489 }
 490
 491 int cv::countNonZero( InputArray _src )
 492 {
 493     Mat src = _src.getMat();
 494     CountNonZeroFunc func = countNonZeroTab[src.depth()];
 495
 496     CV_Assert( src.channels() == 1 && func != 0 );
 497
 498     const Mat* arrays[] = {&src, 0};
 499     uchar* ptrs[1];
 500     NAryMatIterator it(arrays, ptrs);
 501     int total = (int)it.size, nz = 0;
 502
 503     for( size_t i = 0; i < it.nplanes; i++, ++it )
 504         nz += func( ptrs[0], total );
 505
 506     return nz;
 507 }
 508
 509 cv::Scalar cv::mean( InputArray _src, InputArray _mask )
 510 {
 511     Mat src = _src.getMat(), mask = _mask.getMat();
 512     CV_Assert( mask.empty() || mask.type() == CV_8U );
 513
 514     int k, cn = src.channels(), depth = src.depth();
 515     SumFunc func = sumTab[depth];
 516
 517     CV_Assert( cn <= 4 && func != 0 );
 518
 519     const Mat* arrays[] = {&src, &mask, 0};
 520     uchar* ptrs[2];
 521     NAryMatIterator it(arrays, ptrs);
 522     Scalar s;
 523     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
 524     int j, count = 0;
 525     AutoBuffer<int> _buf;
 526     int* buf = (int*)&s[0];
 527     bool blockSum = depth <= CV_16S;
 528     size_t esz = 0, nz0 = 0;
 529
 530     if( blockSum )
 531     {
 532         intSumBlockSize = depth <= CV_8S ? (1 << 23) : (1 << 15);
 533         blockSize = std::min(blockSize, intSumBlockSize);
 534         _buf.allocate(cn);
 535         buf = _buf;
 536
 537         for( k = 0; k < cn; k++ )
 538             buf[k] = 0;
 539         esz = src.elemSize();
 540     }
 541
 542     for( size_t i = 0; i < it.nplanes; i++, ++it )
 543     {
 544         for( j = 0; j < total; j += blockSize )
 545         {
 546             int bsz = std::min(total - j, blockSize);
 547             int nz = func( ptrs[0], ptrs[1], (uchar*)buf, bsz, cn );
 548             count += nz;
 549             nz0 += nz;
 550             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
 551             {
 552                 for( k = 0; k < cn; k++ )
 553                 {
 554                     s[k] += buf[k];
 555                     buf[k] = 0;
 556                 }
 557                 count = 0;
 558             }
 559             ptrs[0] += bsz*esz;
 560             if( ptrs[1] )
 561                 ptrs[1] += bsz;
 562         }
 563     }
 564     return s*(nz0 ? 1./nz0 : 0);
 565 }
 566
 567
 568 void cv::meanStdDev( InputArray _src, OutputArray _mean, OutputArray _sdv, InputArray _mask )
 569 {
 570     Mat src = _src.getMat(), mask = _mask.getMat();
 571     CV_Assert( mask.empty() || mask.type() == CV_8U );
 572
 573     int k, cn = src.channels(), depth = src.depth();
 574     SumSqrFunc func = sumSqrTab[depth];
 575
 576     CV_Assert( func != 0 );
 577
 578     const Mat* arrays[] = {&src, &mask, 0};
 579     uchar* ptrs[2];
 580     NAryMatIterator it(arrays, ptrs);
 581     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
 582     int j, count = 0, nz0 = 0;
 583     AutoBuffer<double> _buf(cn*4);
 584     double *s = (double*)_buf, *sq = s + cn;
 585     int *sbuf = (int*)s, *sqbuf = (int*)sq;
 586     bool blockSum = depth <= CV_16S, blockSqSum = depth <= CV_8S;
 587     size_t esz = 0;
 588
 589     for( k = 0; k < cn; k++ )
 590         s[k] = sq[k] = 0;
 591
 592     if( blockSum )
 593     {
 594         intSumBlockSize = 1 << 15;
 595         blockSize = std::min(blockSize, intSumBlockSize);
 596         sbuf = (int*)(sq + cn);
 597         if( blockSqSum )
 598             sqbuf = sbuf + cn;
 599         for( k = 0; k < cn; k++ )
 600             sbuf[k] = sqbuf[k] = 0;
 601         esz = src.elemSize();
 602     }
 603
 604     for( size_t i = 0; i < it.nplanes; i++, ++it )
 605     {
 606         for( j = 0; j < total; j += blockSize )
 607         {
 608             int bsz = std::min(total - j, blockSize);
 609             int nz = func( ptrs[0], ptrs[1], (uchar*)sbuf, (uchar*)sqbuf, bsz, cn );
 610             count += nz;
 611             nz0 += nz;
 612             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
 613             {
 614                 for( k = 0; k < cn; k++ )
 615                 {
 616                     s[k] += sbuf[k];
 617                     sbuf[k] = 0;
 618                 }
 619                 if( blockSqSum )
 620                 {
 621                     for( k = 0; k < cn; k++ )
 622                     {
 623                         sq[k] += sqbuf[k];
 624                         sqbuf[k] = 0;
 625                     }
 626                 }
 627                 count = 0;
 628             }
 629             ptrs[0] += bsz*esz;
 630             if( ptrs[1] )
 631                 ptrs[1] += bsz;
 632         }
 633     }
 634
 635     double scale = nz0 ? 1./nz0 : 0.;
 636     for( k = 0; k < cn; k++ )
 637     {
 638         s[k] *= scale;
 639         sq[k] = std::sqrt(std::max(sq[k]*scale - s[k]*s[k], 0.));
 640     }
 641
 642     for( j = 0; j < 2; j++ )
 643     {
 644         const double* sptr = j == 0 ? s : sq;
 645         _OutputArray _dst = j == 0 ? _mean : _sdv;
 646         if( !_dst.needed() )
 647             continue;
 648
 649         if( !_dst.fixedSize() )
 650             _dst.create(cn, 1, CV_64F, -1, true);
 651         Mat dst = _dst.getMat();
 652         int dcn = (int)dst.total();
 653         CV_Assert( dst.type() == CV_64F && dst.isContinuous() &&
 654                    (dst.cols == 1 || dst.rows == 1) && dcn >= cn );
 655         double* dptr = dst.ptr<double>();
 656         for( k = 0; k < cn; k++ )
 657             dptr[k] = sptr[k];
 658         for( ; k < dcn; k++ )
 659             dptr[k] = 0;
 660     }
 661 }
 662
 663 /****************************************************************************************\
 664 *                                       minMaxLoc                                        *
 665 \****************************************************************************************/
 666
 667 namespace cv
 668 {
 669
 670 template<typename T, typename WT> static void
 671 minMaxIdx_( const T* src, const uchar* mask, WT* _minVal, WT* _maxVal,
 672             size_t* _minIdx, size_t* _maxIdx, int len, size_t startIdx )
 673 {
 674     WT minVal = *_minVal, maxVal = *_maxVal;
 675     size_t minIdx = *_minIdx, maxIdx = *_maxIdx;
 676
 677     if( !mask )
 678     {
 679         for( int i = 0; i < len; i++ )
 680         {
 681             T val = src[i];
 682             if( val < minVal )
 683             {
 684                 minVal = val;
 685                 minIdx = startIdx + i;
 686             }
 687             if( val > maxVal )
 688             {
 689                 maxVal = val;
 690                 maxIdx = startIdx + i;
 691             }
 692         }
 693     }
 694     else
 695     {
 696         for( int i = 0; i < len; i++ )
 697         {
 698             T val = src[i];
 699             if( mask[i] && val < minVal )
 700             {
 701                 minVal = val;
 702                 minIdx = startIdx + i;
 703             }
 704             if( mask[i] && val > maxVal )
 705             {
 706                 maxVal = val;
 707                 maxIdx = startIdx + i;
 708             }
 709         }
 710     }
 711
 712     *_minIdx = minIdx;
 713     *_maxIdx = maxIdx;
 714     *_minVal = minVal;
 715     *_maxVal = maxVal;
 716 }
 717
 718 static void minMaxIdx_8u(const uchar* src, const uchar* mask, int* minval, int* maxval,
 719                          size_t* minidx, size_t* maxidx, int len, size_t startidx )
 720 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 721
 722 static void minMaxIdx_8s(const schar* src, const uchar* mask, int* minval, int* maxval,
 723                          size_t* minidx, size_t* maxidx, int len, size_t startidx )
 724 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 725
 726 static void minMaxIdx_16u(const ushort* src, const uchar* mask, int* minval, int* maxval,
 727                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 728 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 729
 730 static void minMaxIdx_16s(const short* src, const uchar* mask, int* minval, int* maxval,
 731                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 732 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 733
 734 static void minMaxIdx_32s(const int* src, const uchar* mask, int* minval, int* maxval,
 735                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 736 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 737
 738 static void minMaxIdx_32f(const float* src, const uchar* mask, float* minval, float* maxval,
 739                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 740 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 741
 742 static void minMaxIdx_64f(const double* src, const uchar* mask, double* minval, double* maxval,
 743                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
 744 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
 745
 746 typedef void (*MinMaxIdxFunc)(const uchar*, const uchar*, int*, int*, size_t*, size_t*, int, size_t);
 747
 748 static MinMaxIdxFunc minmaxTab[] =
 749 {
 750     (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_8u), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_8s),
 751     (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_16u), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_16s),
 752     (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_32s),
 753     (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_32f), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_64f),
 754     0
 755 };
 756
 757 static void ofs2idx(const Mat& a, size_t ofs, int* idx)
 758 {
 759     int i, d = a.dims;
 760     if( ofs > 0 )
 761     {
 762         ofs--;
 763         for( i = d-1; i >= 0; i-- )
 764         {
 765             int sz = a.size[i];
 766             idx[i] = (int)(ofs % sz);
 767             ofs /= sz;
 768         }
 769     }
 770     else
 771     {
 772         for( i = d-1; i >= 0; i-- )
 773             idx[i] = -1;
 774     }
 775 }
 776
 777 }
 778
 779 void cv::minMaxIdx(InputArray _src, double* minVal,
 780                    double* maxVal, int* minIdx, int* maxIdx,
 781                    InputArray _mask)
 782 {
 783     Mat src = _src.getMat(), mask = _mask.getMat();
 784     int depth = src.depth(), cn = src.channels();
 785
 786     CV_Assert( (cn == 1 && (mask.empty() || mask.type() == CV_8U)) ||
 787                (cn >= 1 && mask.empty() && !minIdx && !maxIdx) );
 788     MinMaxIdxFunc func = minmaxTab[depth];
 789     CV_Assert( func != 0 );
 790
 791     const Mat* arrays[] = {&src, &mask, 0};
 792     uchar* ptrs[2];
 793     NAryMatIterator it(arrays, ptrs);
 794
 795     size_t minidx = 0, maxidx = 0;
 796     int iminval = INT_MAX, imaxval = INT_MIN;
 797     float fminval = FLT_MAX, fmaxval = -FLT_MAX;
 798     double dminval = DBL_MAX, dmaxval = -DBL_MAX;
 799     size_t startidx = 1;
 800     int *minval = &iminval, *maxval = &imaxval;
 801     int planeSize = (int)it.size*cn;
 802
 803     if( depth == CV_32F )
 804         minval = (int*)&fminval, maxval = (int*)&fmaxval;
 805     else if( depth == CV_64F )
 806         minval = (int*)&dminval, maxval = (int*)&dmaxval;
 807
 808     for( size_t i = 0; i < it.nplanes; i++, ++it, startidx += planeSize )
 809         func( ptrs[0], ptrs[1], minval, maxval, &minidx, &maxidx, planeSize, startidx );
 810
 811     if( minidx == 0 )
 812         dminval = dmaxval = 0;
 813     else if( depth == CV_32F )
 814         dminval = fminval, dmaxval = fmaxval;
 815     else if( depth <= CV_32S )
 816         dminval = iminval, dmaxval = imaxval;
 817
 818     if( minVal )
 819         *minVal = dminval;
 820     if( maxVal )
 821         *maxVal = dmaxval;
 822
 823     if( minIdx )
 824         ofs2idx(src, minidx, minIdx);
 825     if( maxIdx )
 826         ofs2idx(src, maxidx, maxIdx);
 827 }
 828
 829 void cv::minMaxLoc( InputArray _img, double* minVal, double* maxVal,
 830                     Point* minLoc, Point* maxLoc, InputArray mask )
 831 {
 832     Mat img = _img.getMat();
 833     CV_Assert(img.dims <= 2);
 834
 835     minMaxIdx(_img, minVal, maxVal, (int*)minLoc, (int*)maxLoc, mask);
 836     if( minLoc )
 837         std::swap(minLoc->x, minLoc->y);
 838     if( maxLoc )
 839         std::swap(maxLoc->x, maxLoc->y);
 840 }
 841
 842 /****************************************************************************************\
 843 *                                         norm                                           *
 844 \****************************************************************************************/
 845
 846 namespace cv
 847 {
 848
 849 float normL2Sqr_(const float* a, const float* b, int n)
 850 {
 851     int j = 0; float d = 0.f;
 852 #if CV_SSE
 853     if( USE_SSE2 )
 854     {
 855         float CV_DECL_ALIGNED(16) buf[4];
 856         __m128 d0 = _mm_setzero_ps(), d1 = _mm_setzero_ps();
 857
 858         for( ; j <= n - 8; j += 8 )
 859         {
 860             __m128 t0 = _mm_sub_ps(_mm_loadu_ps(a + j), _mm_loadu_ps(b + j));
 861             __m128 t1 = _mm_sub_ps(_mm_loadu_ps(a + j + 4), _mm_loadu_ps(b + j + 4));
 862             d0 = _mm_add_ps(d0, _mm_mul_ps(t0, t0));
 863             d1 = _mm_add_ps(d1, _mm_mul_ps(t1, t1));
 864         }
 865         _mm_store_ps(buf, _mm_add_ps(d0, d1));
 866         d = buf[0] + buf[1] + buf[2] + buf[3];
 867     }
 868     else
 869 #endif
 870     {
 871         for( ; j <= n - 4; j += 4 )
 872         {
 873             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];
 874             d += t0*t0 + t1*t1 + t2*t2 + t3*t3;
 875         }
 876     }
 877
 878     for( ; j < n; j++ )
 879     {
 880         float t = a[j] - b[j];
 881         d += t*t;
 882     }
 883     return d;
 884 }
 885
 886
 887 float normL1_(const float* a, const float* b, int n)
 888 {
 889     int j = 0; float d = 0.f;
 890 #if CV_SSE
 891     if( USE_SSE2 )
 892     {
 893         float CV_DECL_ALIGNED(16) buf[4];
 894         static const int CV_DECL_ALIGNED(16) absbuf[4] = {0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff};
 895         __m128 d0 = _mm_setzero_ps(), d1 = _mm_setzero_ps();
 896         __m128 absmask = _mm_load_ps((const float*)absbuf);
 897
 898         for( ; j <= n - 8; j += 8 )
 899         {
 900             __m128 t0 = _mm_sub_ps(_mm_loadu_ps(a + j), _mm_loadu_ps(b + j));
 901             __m128 t1 = _mm_sub_ps(_mm_loadu_ps(a + j + 4), _mm_loadu_ps(b + j + 4));
 902             d0 = _mm_add_ps(d0, _mm_and_ps(t0, absmask));
 903             d1 = _mm_add_ps(d1, _mm_and_ps(t1, absmask));
 904         }
 905         _mm_store_ps(buf, _mm_add_ps(d0, d1));
 906         d = buf[0] + buf[1] + buf[2] + buf[3];
 907     }
 908     else
 909 #endif
 910     {
 911         for( ; j <= n - 4; j += 4 )
 912         {
 913             d += std::abs(a[j] - b[j]) + std::abs(a[j+1] - b[j+1]) +
 914                     std::abs(a[j+2] - b[j+2]) + std::abs(a[j+3] - b[j+3]);
 915         }
 916     }
 917
 918     for( ; j < n; j++ )
 919         d += std::abs(a[j] - b[j]);
 920     return d;
 921 }
 922
 923 int normL1_(const uchar* a, const uchar* b, int n)
 924 {
 925     int j = 0, d = 0;
 926 #if CV_SSE
 927     if( USE_SSE2 )
 928     {
 929         __m128i d0 = _mm_setzero_si128();
 930
 931         for( ; j <= n - 16; j += 16 )
 932         {
 933             __m128i t0 = _mm_loadu_si128((const __m128i*)(a + j));
 934             __m128i t1 = _mm_loadu_si128((const __m128i*)(b + j));
 935
 936             d0 = _mm_add_epi32(d0, _mm_sad_epu8(t0, t1));
 937         }
 938
 939         for( ; j <= n - 4; j += 4 )
 940         {
 941             __m128i t0 = _mm_cvtsi32_si128(*(const int*)(a + j));
 942             __m128i t1 = _mm_cvtsi32_si128(*(const int*)(b + j));
 943
 944             d0 = _mm_add_epi32(d0, _mm_sad_epu8(t0, t1));
 945         }
 946         d = _mm_cvtsi128_si32(_mm_add_epi32(d0, _mm_unpackhi_epi64(d0, d0)));
 947     }
 948     else
 949 #endif
 950     {
 951         for( ; j <= n - 4; j += 4 )
 952         {
 953             d += std::abs(a[j] - b[j]) + std::abs(a[j+1] - b[j+1]) +
 954                     std::abs(a[j+2] - b[j+2]) + std::abs(a[j+3] - b[j+3]);
 955         }
 956     }
 957     for( ; j < n; j++ )
 958         d += std::abs(a[j] - b[j]);
 959     return d;
 960 }
 961
 962 static const uchar popCountTable[] =
 963 {
 964     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,
 965     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,
 966     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,
 967     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,
 968     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,
 969     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,
 970     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,
 971     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
 972 };
 973
 974 static const uchar popCountTable2[] =
 975 {
 976     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,
 977     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,
 978     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,
 979     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,
 980     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,
 981     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,
 982     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,
 983     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
 984 };
 985
 986 static const uchar popCountTable4[] =
 987 {
 988     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,
 989     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,
 990     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,
 991     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,
 992     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,
 993     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,
 994     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,
 995     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
 996 };
 997
 998 static int normHamming(const uchar* a, int n)
 999 {
1000     int i = 0, result = 0;
1001 #if CV_NEON
1002     uint32x4_t bits = vmovq_n_u32(0);
1003     for (; i <= n - 16; i += 16) {
1004         uint8x16_t A_vec = vld1q_u8 (a + i);
1005         uint8x16_t bitsSet = vcntq_u8 (A_vec);
1006         uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
1007         uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
1008         bits = vaddq_u32(bits, bitSet4);
1009     }
1010     uint64x2_t bitSet2 = vpaddlq_u32 (bits);
1011     result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
1012     result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
1013 #else
1014     for( ; i <= n - 4; i += 4 )
1015             result += popCountTable[a[i]] + popCountTable[a[i+1]] +
1016             popCountTable[a[i+2]] + popCountTable[a[i+3]];
1017 #endif
1018     for( ; i < n; i++ )
1019         result += popCountTable[a[i]];
1020     return result;
1021 }
1022
1023 int normHamming(const uchar* a, const uchar* b, int n)
1024 {
1025     int i = 0, result = 0;
1026 #if CV_NEON
1027     uint32x4_t bits = vmovq_n_u32(0);
1028     for (; i <= n - 16; i += 16) {
1029         uint8x16_t A_vec = vld1q_u8 (a + i);
1030         uint8x16_t B_vec = vld1q_u8 (b + i);
1031         uint8x16_t AxorB = veorq_u8 (A_vec, B_vec);
1032         uint8x16_t bitsSet = vcntq_u8 (AxorB);
1033         uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
1034         uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
1035         bits = vaddq_u32(bits, bitSet4);
1036     }
1037     uint64x2_t bitSet2 = vpaddlq_u32 (bits);
1038     result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
1039     result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
1040 #else
1041     for( ; i <= n - 4; i += 4 )
1042         result += popCountTable[a[i] ^ b[i]] + popCountTable[a[i+1] ^ b[i+1]] +
1043                 popCountTable[a[i+2] ^ b[i+2]] + popCountTable[a[i+3] ^ b[i+3]];
1044 #endif
1045     for( ; i < n; i++ )
1046         result += popCountTable[a[i] ^ b[i]];
1047     return result;
1048 }
1049
1050 static int normHamming(const uchar* a, int n, int cellSize)
1051 {
1052     if( cellSize == 1 )
1053         return normHamming(a, n);
1054     const uchar* tab = 0;
1055     if( cellSize == 2 )
1056         tab = popCountTable2;
1057     else if( cellSize == 4 )
1058         tab = popCountTable4;
1059     else
1060         CV_Error( CV_StsBadSize, "bad cell size (not 1, 2 or 4) in normHamming" );
1061     int i = 0, result = 0;
1062 #if CV_ENABLE_UNROLLED
1063     for( ; i <= n - 4; i += 4 )
1064         result += tab[a[i]] + tab[a[i+1]] + tab[a[i+2]] + tab[a[i+3]];
1065 #endif
1066     for( ; i < n; i++ )
1067         result += tab[a[i]];
1068     return result;
1069 }
1070
1071 int normHamming(const uchar* a, const uchar* b, int n, int cellSize)
1072 {
1073     if( cellSize == 1 )
1074         return normHamming(a, b, n);
1075     const uchar* tab = 0;
1076     if( cellSize == 2 )
1077         tab = popCountTable2;
1078     else if( cellSize == 4 )
1079         tab = popCountTable4;
1080     else
1081         CV_Error( CV_StsBadSize, "bad cell size (not 1, 2 or 4) in normHamming" );
1082     int i = 0, result = 0;
1083     #if CV_ENABLE_UNROLLED
1084     for( ; i <= n - 4; i += 4 )
1085         result += tab[a[i] ^ b[i]] + tab[a[i+1] ^ b[i+1]] +
1086                 tab[a[i+2] ^ b[i+2]] + tab[a[i+3] ^ b[i+3]];
1087     #endif
1088     for( ; i < n; i++ )
1089         result += tab[a[i] ^ b[i]];
1090     return result;
1091 }
1092
1093
1094 template<typename T, typename ST> int
1095 normInf_(const T* src, const uchar* mask, ST* _result, int len, int cn)
1096 {
1097     ST result = *_result;
1098     if( !mask )
1099     {
1100         result = std::max(result, normInf<T, ST>(src, len*cn));
1101     }
1102     else
1103     {
1104         for( int i = 0; i < len; i++, src += cn )
1105             if( mask[i] )
1106             {
1107                 for( int k = 0; k < cn; k++ )
1108                     result = std::max(result, ST(fast_abs(src[k])));
1109             }
1110     }
1111     *_result = result;
1112     return 0;
1113 }
1114
1115 template<typename T, typename ST> int
1116 normL1_(const T* src, const uchar* mask, ST* _result, int len, int cn)
1117 {
1118     ST result = *_result;
1119     if( !mask )
1120     {
1121         result += normL1<T, ST>(src, len*cn);
1122     }
1123     else
1124     {
1125         for( int i = 0; i < len; i++, src += cn )
1126             if( mask[i] )
1127             {
1128                 for( int k = 0; k < cn; k++ )
1129                     result += fast_abs(src[k]);
1130             }
1131     }
1132     *_result = result;
1133     return 0;
1134 }
1135
1136 template<typename T, typename ST> int
1137 normL2_(const T* src, const uchar* mask, ST* _result, int len, int cn)
1138 {
1139     ST result = *_result;
1140     if( !mask )
1141     {
1142         result += normL2Sqr<T, ST>(src, len*cn);
1143     }
1144     else
1145     {
1146         for( int i = 0; i < len; i++, src += cn )
1147             if( mask[i] )
1148             {
1149                 for( int k = 0; k < cn; k++ )
1150                 {
1151                     T v = src[k];
1152                     result += (ST)v*v;
1153                 }
1154             }
1155     }
1156     *_result = result;
1157     return 0;
1158 }
1159
1160 template<typename T, typename ST> int
1161 normDiffInf_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
1162 {
1163     ST result = *_result;
1164     if( !mask )
1165     {
1166         result = std::max(result, normInf<T, ST>(src1, src2, len*cn));
1167     }
1168     else
1169     {
1170         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
1171             if( mask[i] )
1172             {
1173                 for( int k = 0; k < cn; k++ )
1174                     result = std::max(result, (ST)std::abs(src1[k] - src2[k]));
1175             }
1176     }
1177     *_result = result;
1178     return 0;
1179 }
1180
1181 template<typename T, typename ST> int
1182 normDiffL1_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
1183 {
1184     ST result = *_result;
1185     if( !mask )
1186     {
1187         result += normL1<T, ST>(src1, src2, len*cn);
1188     }
1189     else
1190     {
1191         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
1192             if( mask[i] )
1193             {
1194                 for( int k = 0; k < cn; k++ )
1195                     result += std::abs(src1[k] - src2[k]);
1196             }
1197     }
1198     *_result = result;
1199     return 0;
1200 }
1201
1202 template<typename T, typename ST> int
1203 normDiffL2_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
1204 {
1205     ST result = *_result;
1206     if( !mask )
1207     {
1208         result += normL2Sqr<T, ST>(src1, src2, len*cn);
1209     }
1210     else
1211     {
1212         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
1213             if( mask[i] )
1214             {
1215                 for( int k = 0; k < cn; k++ )
1216                 {
1217                     ST v = src1[k] - src2[k];
1218                     result += v*v;
1219                 }
1220             }
1221     }
1222     *_result = result;
1223     return 0;
1224 }
1225
1226
1227 #define CV_DEF_NORM_FUNC(L, suffix, type, ntype) \
1228     static int norm##L##_##suffix(const type* src, const uchar* mask, ntype* r, int len, int cn) \
1229 { return norm##L##_(src, mask, r, len, cn); } \
1230     static int normDiff##L##_##suffix(const type* src1, const type* src2, \
1231     const uchar* mask, ntype* r, int len, int cn) \
1232 { return normDiff##L##_(src1, src2, mask, r, (int)len, cn); }
1233
1234 #define CV_DEF_NORM_ALL(suffix, type, inftype, l1type, l2type) \
1235     CV_DEF_NORM_FUNC(Inf, suffix, type, inftype) \
1236     CV_DEF_NORM_FUNC(L1, suffix, type, l1type) \
1237     CV_DEF_NORM_FUNC(L2, suffix, type, l2type)
1238
1239 CV_DEF_NORM_ALL(8u, uchar, int, int, int)
1240 CV_DEF_NORM_ALL(8s, schar, int, int, int)
1241 CV_DEF_NORM_ALL(16u, ushort, int, int, double)
1242 CV_DEF_NORM_ALL(16s, short, int, int, double)
1243 CV_DEF_NORM_ALL(32s, int, int, double, double)
1244 CV_DEF_NORM_ALL(32f, float, float, double, double)
1245 CV_DEF_NORM_ALL(64f, double, double, double, double)
1246
1247
1248 typedef int (*NormFunc)(const uchar*, const uchar*, uchar*, int, int);
1249 typedef int (*NormDiffFunc)(const uchar*, const uchar*, const uchar*, uchar*, int, int);
1250
1251 static NormFunc normTab[3][8] =
1252 {
1253     {
1254         (NormFunc)GET_OPTIMIZED(normInf_8u), (NormFunc)GET_OPTIMIZED(normInf_8s), (NormFunc)GET_OPTIMIZED(normInf_16u), (NormFunc)GET_OPTIMIZED(normInf_16s),
1255         (NormFunc)GET_OPTIMIZED(normInf_32s), (NormFunc)GET_OPTIMIZED(normInf_32f), (NormFunc)normInf_64f, 0
1256     },
1257     {
1258         (NormFunc)GET_OPTIMIZED(normL1_8u), (NormFunc)GET_OPTIMIZED(normL1_8s), (NormFunc)GET_OPTIMIZED(normL1_16u), (NormFunc)GET_OPTIMIZED(normL1_16s),
1259         (NormFunc)GET_OPTIMIZED(normL1_32s), (NormFunc)GET_OPTIMIZED(normL1_32f), (NormFunc)normL1_64f, 0
1260     },
1261     {
1262         (NormFunc)GET_OPTIMIZED(normL2_8u), (NormFunc)GET_OPTIMIZED(normL2_8s), (NormFunc)GET_OPTIMIZED(normL2_16u), (NormFunc)GET_OPTIMIZED(normL2_16s),
1263         (NormFunc)GET_OPTIMIZED(normL2_32s), (NormFunc)GET_OPTIMIZED(normL2_32f), (NormFunc)normL2_64f, 0
1264     }
1265 };
1266
1267 static NormDiffFunc normDiffTab[3][8] =
1268 {
1269     {
1270         (NormDiffFunc)GET_OPTIMIZED(normDiffInf_8u), (NormDiffFunc)normDiffInf_8s,
1271         (NormDiffFunc)normDiffInf_16u, (NormDiffFunc)normDiffInf_16s,
1272         (NormDiffFunc)normDiffInf_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffInf_32f),
1273         (NormDiffFunc)normDiffInf_64f, 0
1274     },
1275     {
1276         (NormDiffFunc)GET_OPTIMIZED(normDiffL1_8u), (NormDiffFunc)normDiffL1_8s,
1277         (NormDiffFunc)normDiffL1_16u, (NormDiffFunc)normDiffL1_16s,
1278         (NormDiffFunc)normDiffL1_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffL1_32f),
1279         (NormDiffFunc)normDiffL1_64f, 0
1280     },
1281     {
1282         (NormDiffFunc)GET_OPTIMIZED(normDiffL2_8u), (NormDiffFunc)normDiffL2_8s,
1283         (NormDiffFunc)normDiffL2_16u, (NormDiffFunc)normDiffL2_16s,
1284         (NormDiffFunc)normDiffL2_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffL2_32f),
1285         (NormDiffFunc)normDiffL2_64f, 0
1286     }
1287 };
1288
1289 }
1290
1291 double cv::norm( InputArray _src, int normType, InputArray _mask )
1292 {
1293     Mat src = _src.getMat(), mask = _mask.getMat();
1294     int depth = src.depth(), cn = src.channels();
1295
1296     normType &= 7;
1297     CV_Assert( normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 || normType == NORM_L2SQR ||
1298                ((normType == NORM_HAMMING || normType == NORM_HAMMING2) && src.type() == CV_8U) );
1299
1300     if( src.isContinuous() && mask.empty() )
1301     {
1302         size_t len = src.total()*cn;
1303         if( len == (size_t)(int)len )
1304         {
1305             if( depth == CV_32F )
1306             {
1307                 const float* data = src.ptr<float>();
1308
1309                 if( normType == NORM_L2 )
1310                 {
1311                     double result = 0;
1312                     GET_OPTIMIZED(normL2_32f)(data, 0, &result, (int)len, 1);
1313                     return std::sqrt(result);
1314                 }
1315                 if( normType == NORM_L2SQR )
1316                 {
1317                     double result = 0;
1318                     GET_OPTIMIZED(normL2_32f)(data, 0, &result, (int)len, 1);
1319                     return result;
1320                 }
1321                 if( normType == NORM_L1 )
1322                 {
1323                     double result = 0;
1324                     GET_OPTIMIZED(normL1_32f)(data, 0, &result, (int)len, 1);
1325                     return result;
1326                 }
1327                 if( normType == NORM_INF )
1328                 {
1329                     float result = 0;
1330                     GET_OPTIMIZED(normInf_32f)(data, 0, &result, (int)len, 1);
1331                     return result;
1332                 }
1333             }
1334             if( depth == CV_8U )
1335             {
1336                 const uchar* data = src.ptr<uchar>();
1337
1338                 if( normType == NORM_HAMMING )
1339                     return normHamming(data, (int)len);
1340
1341                 if( normType == NORM_HAMMING2 )
1342                     return normHamming(data, (int)len, 2);
1343             }
1344         }
1345     }
1346
1347     CV_Assert( mask.empty() || mask.type() == CV_8U );
1348
1349     if( normType == NORM_HAMMING || normType == NORM_HAMMING2 )
1350     {
1351         if( !mask.empty() )
1352         {
1353             Mat temp;
1354             bitwise_and(src, mask, temp);
1355             return norm(temp, normType);
1356         }
1357         int cellSize = normType == NORM_HAMMING ? 1 : 2;
1358
1359         const Mat* arrays[] = {&src, 0};
1360         uchar* ptrs[1];
1361         NAryMatIterator it(arrays, ptrs);
1362         int total = (int)it.size;
1363         int result = 0;
1364
1365         for( size_t i = 0; i < it.nplanes; i++, ++it )
1366             result += normHamming(ptrs[0], total, cellSize);
1367
1368         return result;
1369     }
1370
1371     NormFunc func = normTab[normType >> 1][depth];
1372     CV_Assert( func != 0 );
1373
1374     const Mat* arrays[] = {&src, &mask, 0};
1375     uchar* ptrs[2];
1376     union
1377     {
1378         double d;
1379         int i;
1380         float f;
1381     }
1382     result;
1383     result.d = 0;
1384     NAryMatIterator it(arrays, ptrs);
1385     int j, total = (int)it.size, blockSize = total, intSumBlockSize = 0, count = 0;
1386     bool blockSum = (normType == NORM_L1 && depth <= CV_16S) ||
1387             ((normType == NORM_L2 || normType == NORM_L2SQR) && depth <= CV_8S);
1388     int isum = 0;
1389     int *ibuf = &result.i;
1390     size_t esz = 0;
1391
1392     if( blockSum )
1393     {
1394         intSumBlockSize = (normType == NORM_L1 && depth <= CV_8S ? (1 << 23) : (1 << 15))/cn;
1395         blockSize = std::min(blockSize, intSumBlockSize);
1396         ibuf = &isum;
1397         esz = src.elemSize();
1398     }
1399
1400     for( size_t i = 0; i < it.nplanes; i++, ++it )
1401     {
1402         for( j = 0; j < total; j += blockSize )
1403         {
1404             int bsz = std::min(total - j, blockSize);
1405             func( ptrs[0], ptrs[1], (uchar*)ibuf, bsz, cn );
1406             count += bsz;
1407             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
1408             {
1409                 result.d += isum;
1410                 isum = 0;
1411                 count = 0;
1412             }
1413             ptrs[0] += bsz*esz;
1414             if( ptrs[1] )
1415                 ptrs[1] += bsz;
1416         }
1417     }
1418
1419     if( normType == NORM_INF )
1420     {
1421         if( depth == CV_64F )
1422             ;
1423         else if( depth == CV_32F )
1424             result.d = result.f;
1425         else
1426             result.d = result.i;
1427     }
1428     else if( normType == NORM_L2 )
1429         result.d = std::sqrt(result.d);
1430
1431     return result.d;
1432 }
1433
1434
1435 double cv::norm( InputArray _src1, InputArray _src2, int normType, InputArray _mask )
1436 {
1437     if( normType & CV_RELATIVE )
1438         return norm(_src1, _src2, normType & ~CV_RELATIVE, _mask)/(norm(_src2, normType, _mask) + DBL_EPSILON);
1439
1440     Mat src1 = _src1.getMat(), src2 = _src2.getMat(), mask = _mask.getMat();
1441     int depth = src1.depth(), cn = src1.channels();
1442
1443     CV_Assert( src1.size == src2.size && src1.type() == src2.type() );
1444
1445     normType &= 7;
1446     CV_Assert( normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 || normType == NORM_L2SQR ||
1447               ((normType == NORM_HAMMING || normType == NORM_HAMMING2) && src1.type() == CV_8U) );
1448
1449     if( src1.isContinuous() && src2.isContinuous() && mask.empty() )
1450     {
1451         size_t len = src1.total()*src1.channels();
1452         if( len == (size_t)(int)len )
1453         {
1454             if( src1.depth() == CV_32F )
1455             {
1456                 const float* data1 = src1.ptr<float>();
1457                 const float* data2 = src2.ptr<float>();
1458
1459                 if( normType == NORM_L2 )
1460                 {
1461                     double result = 0;
1462                     GET_OPTIMIZED(normDiffL2_32f)(data1, data2, 0, &result, (int)len, 1);
1463                     return std::sqrt(result);
1464                 }
1465                 if( normType == NORM_L2SQR )
1466                 {
1467                     double result = 0;
1468                     GET_OPTIMIZED(normDiffL2_32f)(data1, data2, 0, &result, (int)len, 1);
1469                     return result;
1470                 }
1471                 if( normType == NORM_L1 )
1472                 {
1473                     double result = 0;
1474                     GET_OPTIMIZED(normDiffL1_32f)(data1, data2, 0, &result, (int)len, 1);
1475                     return result;
1476                 }
1477                 if( normType == NORM_INF )
1478                 {
1479                     float result = 0;
1480                     GET_OPTIMIZED(normDiffInf_32f)(data1, data2, 0, &result, (int)len, 1);
1481                     return result;
1482                 }
1483             }
1484         }
1485     }
1486
1487     CV_Assert( mask.empty() || mask.type() == CV_8U );
1488
1489     if( normType == NORM_HAMMING || normType == NORM_HAMMING2 )
1490     {
1491         if( !mask.empty() )
1492         {
1493             Mat temp;
1494             bitwise_xor(src1, src2, temp);
1495             bitwise_and(temp, mask, temp);
1496             return norm(temp, normType);
1497         }
1498         int cellSize = normType == NORM_HAMMING ? 1 : 2;
1499
1500         const Mat* arrays[] = {&src1, &src2, 0};
1501         uchar* ptrs[2];
1502         NAryMatIterator it(arrays, ptrs);
1503         int total = (int)it.size;
1504         int result = 0;
1505
1506         for( size_t i = 0; i < it.nplanes; i++, ++it )
1507             result += normHamming(ptrs[0], ptrs[1], total, cellSize);
1508
1509         return result;
1510     }
1511
1512     NormDiffFunc func = normDiffTab[normType >> 1][depth];
1513     CV_Assert( func != 0 );
1514
1515     const Mat* arrays[] = {&src1, &src2, &mask, 0};
1516     uchar* ptrs[3];
1517     union
1518     {
1519         double d;
1520         float f;
1521         int i;
1522         unsigned u;
1523     }
1524     result;
1525     result.d = 0;
1526     NAryMatIterator it(arrays, ptrs);
1527     int j, total = (int)it.size, blockSize = total, intSumBlockSize = 0, count = 0;
1528     bool blockSum = (normType == NORM_L1 && depth <= CV_16S) ||
1529             ((normType == NORM_L2 || normType == NORM_L2SQR) && depth <= CV_8S);
1530     unsigned isum = 0;
1531     unsigned *ibuf = &result.u;
1532     size_t esz = 0;
1533
1534     if( blockSum )
1535     {
1536         intSumBlockSize = normType == NORM_L1 && depth <= CV_8S ? (1 << 23) : (1 << 15);
1537         blockSize = std::min(blockSize, intSumBlockSize);
1538         ibuf = &isum;
1539         esz = src1.elemSize();
1540     }
1541
1542     for( size_t i = 0; i < it.nplanes; i++, ++it )
1543     {
1544         for( j = 0; j < total; j += blockSize )
1545         {
1546             int bsz = std::min(total - j, blockSize);
1547             func( ptrs[0], ptrs[1], ptrs[2], (uchar*)ibuf, bsz, cn );
1548             count += bsz;
1549             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
1550             {
1551                 result.d += isum;
1552                 isum = 0;
1553                 count = 0;
1554             }
1555             ptrs[0] += bsz*esz;
1556             ptrs[1] += bsz*esz;
1557             if( ptrs[2] )
1558                 ptrs[2] += bsz;
1559         }
1560     }
1561
1562     if( normType == NORM_INF )
1563     {
1564         if( depth == CV_64F )
1565             ;
1566         else if( depth == CV_32F )
1567             result.d = result.f;
1568         else
1569             result.d = result.u;
1570     }
1571     else if( normType == NORM_L2 )
1572         result.d = std::sqrt(result.d);
1573
1574     return result.d;
1575 }
1576
1577
1578 ///////////////////////////////////// batch distance ///////////////////////////////////////
1579
1580 namespace cv
1581 {
1582
1583 template<typename _Tp, typename _Rt>
1584 void batchDistL1_(const _Tp* src1, const _Tp* src2, size_t step2,
1585                   int nvecs, int len, _Rt* dist, const uchar* mask)
1586 {
1587     step2 /= sizeof(src2[0]);
1588     if( !mask )
1589     {
1590         for( int i = 0; i < nvecs; i++ )
1591             dist[i] = normL1<_Tp, _Rt>(src1, src2 + step2*i, len);
1592     }
1593     else
1594     {
1595         _Rt val0 = std::numeric_limits<_Rt>::max();
1596         for( int i = 0; i < nvecs; i++ )
1597             dist[i] = mask[i] ? normL1<_Tp, _Rt>(src1, src2 + step2*i, len) : val0;
1598     }
1599 }
1600
1601 template<typename _Tp, typename _Rt>
1602 void batchDistL2Sqr_(const _Tp* src1, const _Tp* src2, size_t step2,
1603                      int nvecs, int len, _Rt* dist, const uchar* mask)
1604 {
1605     step2 /= sizeof(src2[0]);
1606     if( !mask )
1607     {
1608         for( int i = 0; i < nvecs; i++ )
1609             dist[i] = normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len);
1610     }
1611     else
1612     {
1613         _Rt val0 = std::numeric_limits<_Rt>::max();
1614         for( int i = 0; i < nvecs; i++ )
1615             dist[i] = mask[i] ? normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len) : val0;
1616     }
1617 }
1618
1619 template<typename _Tp, typename _Rt>
1620 void batchDistL2_(const _Tp* src1, const _Tp* src2, size_t step2,
1621                   int nvecs, int len, _Rt* dist, const uchar* mask)
1622 {
1623     step2 /= sizeof(src2[0]);
1624     if( !mask )
1625     {
1626         for( int i = 0; i < nvecs; i++ )
1627             dist[i] = std::sqrt(normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len));
1628     }
1629     else
1630     {
1631         _Rt val0 = std::numeric_limits<_Rt>::max();
1632         for( int i = 0; i < nvecs; i++ )
1633             dist[i] = mask[i] ? std::sqrt(normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len)) : val0;
1634     }
1635 }
1636
1637 static void batchDistHamming(const uchar* src1, const uchar* src2, size_t step2,
1638                              int nvecs, int len, int* dist, const uchar* mask)
1639 {
1640     step2 /= sizeof(src2[0]);
1641     if( !mask )
1642     {
1643         for( int i = 0; i < nvecs; i++ )
1644             dist[i] = normHamming(src1, src2 + step2*i, len);
1645     }
1646     else
1647     {
1648         int val0 = INT_MAX;
1649         for( int i = 0; i < nvecs; i++ )
1650             dist[i] = mask[i] ? normHamming(src1, src2 + step2*i, len) : val0;
1651     }
1652 }
1653
1654 static void batchDistHamming2(const uchar* src1, const uchar* src2, size_t step2,
1655                               int nvecs, int len, int* dist, const uchar* mask)
1656 {
1657     step2 /= sizeof(src2[0]);
1658     if( !mask )
1659     {
1660         for( int i = 0; i < nvecs; i++ )
1661             dist[i] = normHamming(src1, src2 + step2*i, len, 2);
1662     }
1663     else
1664     {
1665         int val0 = INT_MAX;
1666         for( int i = 0; i < nvecs; i++ )
1667             dist[i] = mask[i] ? normHamming(src1, src2 + step2*i, len, 2) : val0;
1668     }
1669 }
1670
1671 static void batchDistL1_8u32s(const uchar* src1, const uchar* src2, size_t step2,
1672                                int nvecs, int len, int* dist, const uchar* mask)
1673 {
1674     batchDistL1_<uchar, int>(src1, src2, step2, nvecs, len, dist, mask);
1675 }
1676
1677 static void batchDistL1_8u32f(const uchar* src1, const uchar* src2, size_t step2,
1678                                int nvecs, int len, float* dist, const uchar* mask)
1679 {
1680     batchDistL1_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
1681 }
1682
1683 static void batchDistL2Sqr_8u32s(const uchar* src1, const uchar* src2, size_t step2,
1684                                   int nvecs, int len, int* dist, const uchar* mask)
1685 {
1686     batchDistL2Sqr_<uchar, int>(src1, src2, step2, nvecs, len, dist, mask);
1687 }
1688
1689 static void batchDistL2Sqr_8u32f(const uchar* src1, const uchar* src2, size_t step2,
1690                                   int nvecs, int len, float* dist, const uchar* mask)
1691 {
1692     batchDistL2Sqr_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
1693 }
1694
1695 static void batchDistL2_8u32f(const uchar* src1, const uchar* src2, size_t step2,
1696                                int nvecs, int len, float* dist, const uchar* mask)
1697 {
1698     batchDistL2_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
1699 }
1700
1701 static void batchDistL1_32f(const float* src1, const float* src2, size_t step2,
1702                              int nvecs, int len, float* dist, const uchar* mask)
1703 {
1704     batchDistL1_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
1705 }
1706
1707 static void batchDistL2Sqr_32f(const float* src1, const float* src2, size_t step2,
1708                                 int nvecs, int len, float* dist, const uchar* mask)
1709 {
1710     batchDistL2Sqr_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
1711 }
1712
1713 static void batchDistL2_32f(const float* src1, const float* src2, size_t step2,
1714                              int nvecs, int len, float* dist, const uchar* mask)
1715 {
1716     batchDistL2_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
1717 }
1718
1719 typedef void (*BatchDistFunc)(const uchar* src1, const uchar* src2, size_t step2,
1720                               int nvecs, int len, uchar* dist, const uchar* mask);
1721
1722
1723 struct BatchDistInvoker : public ParallelLoopBody
1724 {
1725     BatchDistInvoker( const Mat& _src1, const Mat& _src2,
1726                       Mat& _dist, Mat& _nidx, int _K,
1727                       const Mat& _mask, int _update,
1728                       BatchDistFunc _func)
1729     {
1730         src1 = &_src1;
1731         src2 = &_src2;
1732         dist = &_dist;
1733         nidx = &_nidx;
1734         K = _K;
1735         mask = &_mask;
1736         update = _update;
1737         func = _func;
1738     }
1739
1740     void operator()(const Range& range) const
1741     {
1742         AutoBuffer<int> buf(src2->rows);
1743         int* bufptr = buf;
1744
1745         for( int i = range.start; i < range.end; i++ )
1746         {
1747             func(src1->ptr(i), src2->ptr(), src2->step, src2->rows, src2->cols,
1748                  K > 0 ? (uchar*)bufptr : dist->ptr(i), mask->data ? mask->ptr(i) : 0);
1749
1750             if( K > 0 )
1751             {
1752                 int* nidxptr = nidx->ptr<int>(i);
1753                 // since positive float's can be compared just like int's,
1754                 // we handle both CV_32S and CV_32F cases with a single branch
1755                 int* distptr = (int*)dist->ptr(i);
1756
1757                 int j, k;
1758
1759                 for( j = 0; j < src2->rows; j++ )
1760                 {
1761                     int d = bufptr[j];
1762                     if( d < distptr[K-1] )
1763                     {
1764                         for( k = K-2; k >= 0 && distptr[k] > d; k-- )
1765                         {
1766                             nidxptr[k+1] = nidxptr[k];
1767                             distptr[k+1] = distptr[k];
1768                         }
1769                         nidxptr[k+1] = j + update;
1770                         distptr[k+1] = d;
1771                     }
1772                 }
1773             }
1774         }
1775     }
1776
1777     const Mat *src1;
1778     const Mat *src2;
1779     Mat *dist;
1780     Mat *nidx;
1781     const Mat *mask;
1782     int K;
1783     int update;
1784     BatchDistFunc func;
1785 };
1786
1787 }
1788
1789 void cv::batchDistance( InputArray _src1, InputArray _src2,
1790                         OutputArray _dist, int dtype, OutputArray _nidx,
1791                         int normType, int K, InputArray _mask,
1792                         int update, bool crosscheck )
1793 {
1794     Mat src1 = _src1.getMat(), src2 = _src2.getMat(), mask = _mask.getMat();
1795     int type = src1.type();
1796     CV_Assert( type == src2.type() && src1.cols == src2.cols &&
1797                (type == CV_32F || type == CV_8U));
1798     CV_Assert( _nidx.needed() == (K > 0) );
1799
1800     if( dtype == -1 )
1801     {
1802         dtype = normType == NORM_HAMMING || normType == NORM_HAMMING2 ? CV_32S : CV_32F;
1803     }
1804     CV_Assert( (type == CV_8U && dtype == CV_32S) || dtype == CV_32F);
1805
1806     K = std::min(K, src2.rows);
1807
1808     _dist.create(src1.rows, (K > 0 ? K : src2.rows), dtype);
1809     Mat dist = _dist.getMat(), nidx;
1810     if( _nidx.needed() )
1811     {
1812         _nidx.create(dist.size(), CV_32S);
1813         nidx = _nidx.getMat();
1814     }
1815
1816     if( update == 0 && K > 0 )
1817     {
1818         dist = Scalar::all(dtype == CV_32S ? (double)INT_MAX : (double)FLT_MAX);
1819         nidx = Scalar::all(-1);
1820     }
1821
1822     if( crosscheck )
1823     {
1824         CV_Assert( K == 1 && update == 0 && mask.empty() );
1825         Mat tdist, tidx;
1826         batchDistance(src2, src1, tdist, dtype, tidx, normType, K, mask, 0, false);
1827
1828         // if an idx-th element from src1 appeared to be the nearest to i-th element of src2,
1829         // we update the minimum mutual distance between idx-th element of src1 and the whole src2 set.
1830         // As a result, if nidx[idx] = i*, it means that idx-th element of src1 is the nearest
1831         // to i*-th element of src2 and i*-th element of src2 is the closest to idx-th element of src1.
1832         // If nidx[idx] = -1, it means that there is no such ideal couple for it in src2.
1833         // This O(N) procedure is called cross-check and it helps to eliminate some false matches.
1834         if( dtype == CV_32S )
1835         {
1836             for( int i = 0; i < tdist.rows; i++ )
1837             {
1838                 int idx = tidx.at<int>(i);
1839                 int d = tdist.at<int>(i), d0 = dist.at<int>(idx);
1840                 if( d < d0 )
1841                 {
1842                     dist.at<int>(idx) = d;
1843                     nidx.at<int>(idx) = i + update;
1844                 }
1845             }
1846         }
1847         else
1848         {
1849             for( int i = 0; i < tdist.rows; i++ )
1850             {
1851                 int idx = tidx.at<int>(i);
1852                 float d = tdist.at<float>(i), d0 = dist.at<float>(idx);
1853                 if( d < d0 )
1854                 {
1855                     dist.at<float>(idx) = d;
1856                     nidx.at<int>(idx) = i + update;
1857                 }
1858             }
1859         }
1860         return;
1861     }
1862
1863     BatchDistFunc func = 0;
1864     if( type == CV_8U )
1865     {
1866         if( normType == NORM_L1 && dtype == CV_32S )
1867             func = (BatchDistFunc)batchDistL1_8u32s;
1868         else if( normType == NORM_L1 && dtype == CV_32F )
1869             func = (BatchDistFunc)batchDistL1_8u32f;
1870         else if( normType == NORM_L2SQR && dtype == CV_32S )
1871             func = (BatchDistFunc)batchDistL2Sqr_8u32s;
1872         else if( normType == NORM_L2SQR && dtype == CV_32F )
1873             func = (BatchDistFunc)batchDistL2Sqr_8u32f;
1874         else if( normType == NORM_L2 && dtype == CV_32F )
1875             func = (BatchDistFunc)batchDistL2_8u32f;
1876         else if( normType == NORM_HAMMING && dtype == CV_32S )
1877             func = (BatchDistFunc)batchDistHamming;
1878         else if( normType == NORM_HAMMING2 && dtype == CV_32S )
1879             func = (BatchDistFunc)batchDistHamming2;
1880     }
1881     else if( type == CV_32F && dtype == CV_32F )
1882     {
1883         if( normType == NORM_L1 )
1884             func = (BatchDistFunc)batchDistL1_32f;
1885         else if( normType == NORM_L2SQR )
1886             func = (BatchDistFunc)batchDistL2Sqr_32f;
1887         else if( normType == NORM_L2 )
1888             func = (BatchDistFunc)batchDistL2_32f;
1889     }
1890
1891     if( func == 0 )
1892         CV_Error_(CV_StsUnsupportedFormat,
1893                   ("The combination of type=%d, dtype=%d and normType=%d is not supported",
1894                    type, dtype, normType));
1895
1896     parallel_for_(Range(0, src1.rows),
1897                   BatchDistInvoker(src1, src2, dist, nidx, K, mask, update, func));
1898 }
1899
1900
1901 void cv::findNonZero( InputArray _src, OutputArray _idx )
1902 {
1903     Mat src = _src.getMat();
1904     CV_Assert( src.type() == CV_8UC1 );
1905     int n = countNonZero(src);
1906     if( _idx.kind() == _InputArray::MAT && !_idx.getMatRef().isContinuous() )
1907         _idx.release();
1908     _idx.create(n, 1, CV_32SC2);
1909     Mat idx = _idx.getMat();
1910     CV_Assert(idx.isContinuous());
1911     Point* idx_ptr = (Point*)idx.data;
1912
1913     for( int i = 0; i < src.rows; i++ )
1914     {
1915         const uchar* bin_ptr = src.ptr(i);
1916         for( int j = 0; j < src.cols; j++ )
1917             if( bin_ptr[j] )
1918                 *idx_ptr++ = Point(j, i);
1919     }
1920 }
1921
1922
1923 CV_IMPL CvScalar cvSum( const CvArr* srcarr )
1924 {
1925     cv::Scalar sum = cv::sum(cv::cvarrToMat(srcarr, false, true, 1));
1926     if( CV_IS_IMAGE(srcarr) )
1927     {
1928         int coi = cvGetImageCOI((IplImage*)srcarr);
1929         if( coi )
1930         {
1931             CV_Assert( 0 < coi && coi <= 4 );
1932             sum = cv::Scalar(sum[coi-1]);
1933         }
1934     }
1935     return sum;
1936 }
1937
1938 CV_IMPL int cvCountNonZero( const CvArr* imgarr )
1939 {
1940     cv::Mat img = cv::cvarrToMat(imgarr, false, true, 1);
1941     if( img.channels() > 1 )
1942         cv::extractImageCOI(imgarr, img);
1943     return countNonZero(img);
1944 }
1945
1946
1947 CV_IMPL  CvScalar
1948 cvAvg( const void* imgarr, const void* maskarr )
1949 {
1950     cv::Mat img = cv::cvarrToMat(imgarr, false, true, 1);
1951     cv::Scalar mean = !maskarr ? cv::mean(img) : cv::mean(img, cv::cvarrToMat(maskarr));
1952     if( CV_IS_IMAGE(imgarr) )
1953     {
1954         int coi = cvGetImageCOI((IplImage*)imgarr);
1955         if( coi )
1956         {
1957             CV_Assert( 0 < coi && coi <= 4 );
1958             mean = cv::Scalar(mean[coi-1]);
1959         }
1960     }
1961     return mean;
1962 }
1963
1964
1965 CV_IMPL  void
1966 cvAvgSdv( const CvArr* imgarr, CvScalar* _mean, CvScalar* _sdv, const void* maskarr )
1967 {
1968     cv::Scalar mean, sdv;
1969
1970     cv::Mat mask;
1971     if( maskarr )
1972         mask = cv::cvarrToMat(maskarr);
1973
1974     cv::meanStdDev(cv::cvarrToMat(imgarr, false, true, 1), mean, sdv, mask );
1975
1976     if( CV_IS_IMAGE(imgarr) )
1977     {
1978         int coi = cvGetImageCOI((IplImage*)imgarr);
1979         if( coi )
1980         {
1981             CV_Assert( 0 < coi && coi <= 4 );
1982             mean = cv::Scalar(mean[coi-1]);
1983             sdv = cv::Scalar(sdv[coi-1]);
1984         }
1985     }
1986
1987     if( _mean )
1988         *(cv::Scalar*)_mean = mean;
1989     if( _sdv )
1990         *(cv::Scalar*)_sdv = sdv;
1991 }
1992
1993
1994 CV_IMPL void
1995 cvMinMaxLoc( const void* imgarr, double* _minVal, double* _maxVal,
1996              CvPoint* _minLoc, CvPoint* _maxLoc, const void* maskarr )
1997 {
1998     cv::Mat mask, img = cv::cvarrToMat(imgarr, false, true, 1);
1999     if( maskarr )
2000         mask = cv::cvarrToMat(maskarr);
2001     if( img.channels() > 1 )
2002         cv::extractImageCOI(imgarr, img);
2003
2004     cv::minMaxLoc( img, _minVal, _maxVal,
2005                    (cv::Point*)_minLoc, (cv::Point*)_maxLoc, mask );
2006 }
2007
2008
2009 CV_IMPL  double
2010 cvNorm( const void* imgA, const void* imgB, int normType, const void* maskarr )
2011 {
2012     cv::Mat a, mask;
2013     if( !imgA )
2014     {
2015         imgA = imgB;
2016         imgB = 0;
2017     }
2018
2019     a = cv::cvarrToMat(imgA, false, true, 1);
2020     if( maskarr )
2021         mask = cv::cvarrToMat(maskarr);
2022
2023     if( a.channels() > 1 && CV_IS_IMAGE(imgA) && cvGetImageCOI((const IplImage*)imgA) > 0 )
2024         cv::extractImageCOI(imgA, a);
2025
2026     if( !imgB )
2027         return !maskarr ? cv::norm(a, normType) : cv::norm(a, normType, mask);
2028
2029     cv::Mat b = cv::cvarrToMat(imgB, false, true, 1);
2030     if( b.channels() > 1 && CV_IS_IMAGE(imgB) && cvGetImageCOI((const IplImage*)imgB) > 0 )
2031         cv::extractImageCOI(imgB, b);
2032
2033     return !maskarr ? cv::norm(a, b, normType) : cv::norm(a, b, normType, mask);
2034 }