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 #include <limits>
  46
  47 #include "opencl_kernels.hpp"
  48
  49 namespace cv
  50 {
  51
  52 template<typename T> static inline Scalar rawToScalar(const T& v)
  53 {
  54     Scalar s;
  55     typedef typename DataType<T>::channel_type T1;
  56     int i, n = DataType<T>::channels;
  57     for( i = 0; i < n; i++ )
  58         s.val[i] = ((T1*)&v)[i];
  59     return s;
  60 }
  61
  62 /****************************************************************************************\
  63 *                                        sum                                             *
  64 \****************************************************************************************/
  65
  66 template<typename T, typename ST>
  67 static int sum_(const T* src0, const uchar* mask, ST* dst, int len, int cn )
  68 {
  69     const T* src = src0;
  70     if( !mask )
  71     {
  72         int i=0;
  73         int k = cn % 4;
  74         if( k == 1 )
  75         {
  76             ST s0 = dst[0];
  77
  78             #if CV_ENABLE_UNROLLED
  79             for(; i <= len - 4; i += 4, src += cn*4 )
  80                 s0 += src[0] + src[cn] + src[cn*2] + src[cn*3];
  81             #endif
  82             for( ; i < len; i++, src += cn )
  83                 s0 += src[0];
  84             dst[0] = s0;
  85         }
  86         else if( k == 2 )
  87         {
  88             ST s0 = dst[0], s1 = dst[1];
  89             for( i = 0; i < len; i++, src += cn )
  90             {
  91                 s0 += src[0];
  92                 s1 += src[1];
  93             }
  94             dst[0] = s0;
  95             dst[1] = s1;
  96         }
  97         else if( k == 3 )
  98         {
  99             ST s0 = dst[0], s1 = dst[1], s2 = dst[2];
 100             for( i = 0; i < len; i++, src += cn )
 101             {
 102                 s0 += src[0];
 103                 s1 += src[1];
 104                 s2 += src[2];
 105             }
 106             dst[0] = s0;
 107             dst[1] = s1;
 108             dst[2] = s2;
 109         }
 110
 111         for( ; k < cn; k += 4 )
 112         {
 113             src = src0 + k;
 114             ST s0 = dst[k], s1 = dst[k+1], s2 = dst[k+2], s3 = dst[k+3];
 115             for( i = 0; i < len; i++, src += cn )
 116             {
 117                 s0 += src[0]; s1 += src[1];
 118                 s2 += src[2]; s3 += src[3];
 119             }
 120             dst[k] = s0;
 121             dst[k+1] = s1;
 122             dst[k+2] = s2;
 123             dst[k+3] = s3;
 124         }
 125         return len;
 126     }
 127
 128     int i, nzm = 0;
 129     if( cn == 1 )
 130     {
 131         ST s = dst[0];
 132         for( i = 0; i < len; i++ )
 133             if( mask[i] )
 134             {
 135                 s += src[i];
 136                 nzm++;
 137             }
 138         dst[0] = s;
 139     }
 140     else if( cn == 3 )
 141     {
 142         ST s0 = dst[0], s1 = dst[1], s2 = dst[2];
 143         for( i = 0; i < len; i++, src += 3 )
 144             if( mask[i] )
 145             {
 146                 s0 += src[0];
 147                 s1 += src[1];
 148                 s2 += src[2];
 149                 nzm++;
 150             }
 151         dst[0] = s0;
 152         dst[1] = s1;
 153         dst[2] = s2;
 154     }
 155     else
 156     {
 157         for( i = 0; i < len; i++, src += cn )
 158             if( mask[i] )
 159             {
 160                 int k = 0;
 161                 #if CV_ENABLE_UNROLLED
 162                 for( ; k <= cn - 4; k += 4 )
 163                 {
 164                     ST s0, s1;
 165                     s0 = dst[k] + src[k];
 166                     s1 = dst[k+1] + src[k+1];
 167                     dst[k] = s0; dst[k+1] = s1;
 168                     s0 = dst[k+2] + src[k+2];
 169                     s1 = dst[k+3] + src[k+3];
 170                     dst[k+2] = s0; dst[k+3] = s1;
 171                 }
 172                 #endif
 173                 for( ; k < cn; k++ )
 174                     dst[k] += src[k];
 175                 nzm++;
 176             }
 177     }
 178     return nzm;
 179 }
 180
 181
 182 static int sum8u( const uchar* src, const uchar* mask, int* dst, int len, int cn )
 183 { return sum_(src, mask, dst, len, cn); }
 184
 185 static int sum8s( const schar* src, const uchar* mask, int* dst, int len, int cn )
 186 { return sum_(src, mask, dst, len, cn); }
 187
 188 static int sum16u( const ushort* src, const uchar* mask, int* dst, int len, int cn )
 189 { return sum_(src, mask, dst, len, cn); }
 190
 191 static int sum16s( const short* src, const uchar* mask, int* dst, int len, int cn )
 192 { return sum_(src, mask, dst, len, cn); }
 193
 194 static int sum32s( const int* src, const uchar* mask, double* dst, int len, int cn )
 195 { return sum_(src, mask, dst, len, cn); }
 196
 197 static int sum32f( const float* src, const uchar* mask, double* dst, int len, int cn )
 198 { return sum_(src, mask, dst, len, cn); }
 199
 200 static int sum64f( const double* src, const uchar* mask, double* dst, int len, int cn )
 201 { return sum_(src, mask, dst, len, cn); }
 202
 203 typedef int (*SumFunc)(const uchar*, const uchar* mask, uchar*, int, int);
 204
 205 static SumFunc getSumFunc(int depth)
 206 {
 207     static SumFunc sumTab[] =
 208     {
 209         (SumFunc)GET_OPTIMIZED(sum8u), (SumFunc)sum8s,
 210         (SumFunc)sum16u, (SumFunc)sum16s,
 211         (SumFunc)sum32s,
 212         (SumFunc)GET_OPTIMIZED(sum32f), (SumFunc)sum64f,
 213         0
 214     };
 215
 216     return sumTab[depth];
 217 }
 218
 219 template<typename T>
 220 static int countNonZero_(const T* src, int len )
 221 {
 222     int i=0, nz = 0;
 223     #if CV_ENABLE_UNROLLED
 224     for(; i <= len - 4; i += 4 )
 225         nz += (src[i] != 0) + (src[i+1] != 0) + (src[i+2] != 0) + (src[i+3] != 0);
 226     #endif
 227     for( ; i < len; i++ )
 228         nz += src[i] != 0;
 229     return nz;
 230 }
 231
 232 static int countNonZero8u( const uchar* src, int len )
 233 {
 234     int i=0, nz = 0;
 235 #if CV_SSE2
 236     if(USE_SSE2)//5x-6x
 237     {
 238         __m128i pattern = _mm_setzero_si128 ();
 239         static uchar tab[256];
 240         static volatile bool initialized = false;
 241         if( !initialized )
 242         {
 243             // we compute inverse popcount table,
 244             // since we pass (img[x] == 0) mask as index in the table.
 245             for( int j = 0; j < 256; j++ )
 246             {
 247                 int val = 0;
 248                 for( int mask = 1; mask < 256; mask += mask )
 249                     val += (j & mask) == 0;
 250                 tab[j] = (uchar)val;
 251             }
 252             initialized = true;
 253         }
 254
 255         for (; i<=len-16; i+=16)
 256         {
 257             __m128i r0 = _mm_loadu_si128((const __m128i*)(src+i));
 258             int val = _mm_movemask_epi8(_mm_cmpeq_epi8(r0, pattern));
 259             nz += tab[val & 255] + tab[val >> 8];
 260         }
 261     }
 262 #endif
 263     for( ; i < len; i++ )
 264         nz += src[i] != 0;
 265     return nz;
 266 }
 267
 268 static int countNonZero16u( const ushort* src, int len )
 269 { return countNonZero_(src, len); }
 270
 271 static int countNonZero32s( const int* src, int len )
 272 { return countNonZero_(src, len); }
 273
 274 static int countNonZero32f( const float* src, int len )
 275 { return countNonZero_(src, len); }
 276
 277 static int countNonZero64f( const double* src, int len )
 278 { return countNonZero_(src, len); }
 279
 280 typedef int (*CountNonZeroFunc)(const uchar*, int);
 281
 282 static CountNonZeroFunc getCountNonZeroTab(int depth)
 283 {
 284     static CountNonZeroFunc countNonZeroTab[] =
 285     {
 286         (CountNonZeroFunc)GET_OPTIMIZED(countNonZero8u), (CountNonZeroFunc)GET_OPTIMIZED(countNonZero8u),
 287         (CountNonZeroFunc)GET_OPTIMIZED(countNonZero16u), (CountNonZeroFunc)GET_OPTIMIZED(countNonZero16u),
 288         (CountNonZeroFunc)GET_OPTIMIZED(countNonZero32s), (CountNonZeroFunc)GET_OPTIMIZED(countNonZero32f),
 289         (CountNonZeroFunc)GET_OPTIMIZED(countNonZero64f), 0
 290     };
 291
 292     return countNonZeroTab[depth];
 293 }
 294
 295 template<typename T, typename ST, typename SQT>
 296 static int sumsqr_(const T* src0, const uchar* mask, ST* sum, SQT* sqsum, int len, int cn )
 297 {
 298     const T* src = src0;
 299
 300     if( !mask )
 301     {
 302         int i;
 303         int k = cn % 4;
 304
 305         if( k == 1 )
 306         {
 307             ST s0 = sum[0];
 308             SQT sq0 = sqsum[0];
 309             for( i = 0; i < len; i++, src += cn )
 310             {
 311                 T v = src[0];
 312                 s0 += v; sq0 += (SQT)v*v;
 313             }
 314             sum[0] = s0;
 315             sqsum[0] = sq0;
 316         }
 317         else if( k == 2 )
 318         {
 319             ST s0 = sum[0], s1 = sum[1];
 320             SQT sq0 = sqsum[0], sq1 = sqsum[1];
 321             for( i = 0; i < len; i++, src += cn )
 322             {
 323                 T v0 = src[0], v1 = src[1];
 324                 s0 += v0; sq0 += (SQT)v0*v0;
 325                 s1 += v1; sq1 += (SQT)v1*v1;
 326             }
 327             sum[0] = s0; sum[1] = s1;
 328             sqsum[0] = sq0; sqsum[1] = sq1;
 329         }
 330         else if( k == 3 )
 331         {
 332             ST s0 = sum[0], s1 = sum[1], s2 = sum[2];
 333             SQT sq0 = sqsum[0], sq1 = sqsum[1], sq2 = sqsum[2];
 334             for( i = 0; i < len; i++, src += cn )
 335             {
 336                 T v0 = src[0], v1 = src[1], v2 = src[2];
 337                 s0 += v0; sq0 += (SQT)v0*v0;
 338                 s1 += v1; sq1 += (SQT)v1*v1;
 339                 s2 += v2; sq2 += (SQT)v2*v2;
 340             }
 341             sum[0] = s0; sum[1] = s1; sum[2] = s2;
 342             sqsum[0] = sq0; sqsum[1] = sq1; sqsum[2] = sq2;
 343         }
 344
 345         for( ; k < cn; k += 4 )
 346         {
 347             src = src0 + k;
 348             ST s0 = sum[k], s1 = sum[k+1], s2 = sum[k+2], s3 = sum[k+3];
 349             SQT sq0 = sqsum[k], sq1 = sqsum[k+1], sq2 = sqsum[k+2], sq3 = sqsum[k+3];
 350             for( i = 0; i < len; i++, src += cn )
 351             {
 352                 T v0, v1;
 353                 v0 = src[0], v1 = src[1];
 354                 s0 += v0; sq0 += (SQT)v0*v0;
 355                 s1 += v1; sq1 += (SQT)v1*v1;
 356                 v0 = src[2], v1 = src[3];
 357                 s2 += v0; sq2 += (SQT)v0*v0;
 358                 s3 += v1; sq3 += (SQT)v1*v1;
 359             }
 360             sum[k] = s0; sum[k+1] = s1;
 361             sum[k+2] = s2; sum[k+3] = s3;
 362             sqsum[k] = sq0; sqsum[k+1] = sq1;
 363             sqsum[k+2] = sq2; sqsum[k+3] = sq3;
 364         }
 365         return len;
 366     }
 367
 368     int i, nzm = 0;
 369
 370     if( cn == 1 )
 371     {
 372         ST s0 = sum[0];
 373         SQT sq0 = sqsum[0];
 374         for( i = 0; i < len; i++ )
 375             if( mask[i] )
 376             {
 377                 T v = src[i];
 378                 s0 += v; sq0 += (SQT)v*v;
 379                 nzm++;
 380             }
 381         sum[0] = s0;
 382         sqsum[0] = sq0;
 383     }
 384     else if( cn == 3 )
 385     {
 386         ST s0 = sum[0], s1 = sum[1], s2 = sum[2];
 387         SQT sq0 = sqsum[0], sq1 = sqsum[1], sq2 = sqsum[2];
 388         for( i = 0; i < len; i++, src += 3 )
 389             if( mask[i] )
 390             {
 391                 T v0 = src[0], v1 = src[1], v2 = src[2];
 392                 s0 += v0; sq0 += (SQT)v0*v0;
 393                 s1 += v1; sq1 += (SQT)v1*v1;
 394                 s2 += v2; sq2 += (SQT)v2*v2;
 395                 nzm++;
 396             }
 397         sum[0] = s0; sum[1] = s1; sum[2] = s2;
 398         sqsum[0] = sq0; sqsum[1] = sq1; sqsum[2] = sq2;
 399     }
 400     else
 401     {
 402         for( i = 0; i < len; i++, src += cn )
 403             if( mask[i] )
 404             {
 405                 for( int k = 0; k < cn; k++ )
 406                 {
 407                     T v = src[k];
 408                     ST s = sum[k] + v;
 409                     SQT sq = sqsum[k] + (SQT)v*v;
 410                     sum[k] = s; sqsum[k] = sq;
 411                 }
 412                 nzm++;
 413             }
 414     }
 415     return nzm;
 416 }
 417
 418
 419 static int sqsum8u( const uchar* src, const uchar* mask, int* sum, int* sqsum, int len, int cn )
 420 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 421
 422 static int sqsum8s( const schar* src, const uchar* mask, int* sum, int* sqsum, int len, int cn )
 423 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 424
 425 static int sqsum16u( const ushort* src, const uchar* mask, int* sum, double* sqsum, int len, int cn )
 426 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 427
 428 static int sqsum16s( const short* src, const uchar* mask, int* sum, double* sqsum, int len, int cn )
 429 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 430
 431 static int sqsum32s( const int* src, const uchar* mask, double* sum, double* sqsum, int len, int cn )
 432 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 433
 434 static int sqsum32f( const float* src, const uchar* mask, double* sum, double* sqsum, int len, int cn )
 435 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 436
 437 static int sqsum64f( const double* src, const uchar* mask, double* sum, double* sqsum, int len, int cn )
 438 { return sumsqr_(src, mask, sum, sqsum, len, cn); }
 439
 440 typedef int (*SumSqrFunc)(const uchar*, const uchar* mask, uchar*, uchar*, int, int);
 441
 442 static SumSqrFunc getSumSqrTab(int depth)
 443 {
 444     static SumSqrFunc sumSqrTab[] =
 445     {
 446         (SumSqrFunc)GET_OPTIMIZED(sqsum8u), (SumSqrFunc)sqsum8s, (SumSqrFunc)sqsum16u, (SumSqrFunc)sqsum16s,
 447         (SumSqrFunc)sqsum32s, (SumSqrFunc)GET_OPTIMIZED(sqsum32f), (SumSqrFunc)sqsum64f, 0
 448     };
 449
 450     return sumSqrTab[depth];
 451 }
 452
 453 #ifdef HAVE_OPENCL
 454
 455 template <typename T> Scalar ocl_part_sum(Mat m)
 456 {
 457     CV_Assert(m.rows == 1);
 458
 459     Scalar s = Scalar::all(0);
 460     int cn = m.channels();
 461     const T * const ptr = m.ptr<T>(0);
 462
 463     for (int x = 0, w = m.cols * cn; x < w; )
 464         for (int c = 0; c < cn; ++c, ++x)
 465             s[c] += ptr[x];
 466
 467     return s;
 468 }
 469
 470 enum { OCL_OP_SUM = 0, OCL_OP_SUM_ABS =  1, OCL_OP_SUM_SQR = 2 };
 471
 472 static bool ocl_sum( InputArray _src, Scalar & res, int sum_op, InputArray _mask = noArray(),
 473                      InputArray _src2 = noArray(), bool calc2 = false, const Scalar & res2 = Scalar() )
 474 {
 475     CV_Assert(sum_op == OCL_OP_SUM || sum_op == OCL_OP_SUM_ABS || sum_op == OCL_OP_SUM_SQR);
 476
 477     const ocl::Device & dev = ocl::Device::getDefault();
 478     bool doubleSupport = dev.doubleFPConfig() > 0,
 479         haveMask = _mask.kind() != _InputArray::NONE,
 480         haveSrc2 = _src2.kind() != _InputArray::NONE;
 481     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type),
 482             kercn = cn == 1 && !haveMask ? ocl::predictOptimalVectorWidth(_src) : 1,
 483             mcn = std::max(cn, kercn);
 484     CV_Assert(!haveSrc2 || _src2.type() == type);
 485
 486     if ( (!doubleSupport && depth == CV_64F) || cn > 4 )
 487         return false;
 488
 489     int ngroups = dev.maxComputeUnits(), dbsize = ngroups * (calc2 ? 2 : 1);
 490     size_t wgs = dev.maxWorkGroupSize();
 491
 492     int ddepth = std::max(sum_op == OCL_OP_SUM_SQR ? CV_32F : CV_32S, depth),
 493             dtype = CV_MAKE_TYPE(ddepth, cn);
 494     CV_Assert(!haveMask || _mask.type() == CV_8UC1);
 495
 496     int wgs2_aligned = 1;
 497     while (wgs2_aligned < (int)wgs)
 498         wgs2_aligned <<= 1;
 499     wgs2_aligned >>= 1;
 500
 501     static const char * const opMap[3] = { "OP_SUM", "OP_SUM_ABS", "OP_SUM_SQR" };
 502     char cvt[40];
 503     String opts = format("-D srcT=%s -D srcT1=%s -D dstT=%s -D dstTK=%s -D dstT1=%s -D ddepth=%d -D cn=%d"
 504                          " -D convertToDT=%s -D %s -D WGS=%d -D WGS2_ALIGNED=%d%s%s%s%s -D kercn=%d%s%s%s",
 505                          ocl::typeToStr(CV_MAKE_TYPE(depth, mcn)), ocl::typeToStr(depth),
 506                          ocl::typeToStr(dtype), ocl::typeToStr(CV_MAKE_TYPE(ddepth, mcn)),
 507                          ocl::typeToStr(ddepth), ddepth, cn,
 508                          ocl::convertTypeStr(depth, ddepth, mcn, cvt),
 509                          opMap[sum_op], (int)wgs, wgs2_aligned,
 510                          doubleSupport ? " -D DOUBLE_SUPPORT" : "",
 511                          haveMask ? " -D HAVE_MASK" : "",
 512                          _src.isContinuous() ? " -D HAVE_SRC_CONT" : "",
 513                          haveMask && _mask.isContinuous() ? " -D HAVE_MASK_CONT" : "", kercn,
 514                          haveSrc2 ? " -D HAVE_SRC2" : "", calc2 ? " -D OP_CALC2" : "",
 515                          haveSrc2 && _src2.isContinuous() ? " -D HAVE_SRC2_CONT" : "");
 516
 517     ocl::Kernel k("reduce", ocl::core::reduce_oclsrc, opts);
 518     if (k.empty())
 519         return false;
 520
 521     UMat src = _src.getUMat(), src2 = _src2.getUMat(),
 522         db(1, dbsize, dtype), mask = _mask.getUMat();
 523
 524     ocl::KernelArg srcarg = ocl::KernelArg::ReadOnlyNoSize(src),
 525             dbarg = ocl::KernelArg::PtrWriteOnly(db),
 526             maskarg = ocl::KernelArg::ReadOnlyNoSize(mask),
 527             src2arg = ocl::KernelArg::ReadOnlyNoSize(src2);
 528
 529     if (haveMask)
 530     {
 531         if (haveSrc2)
 532             k.args(srcarg, src.cols, (int)src.total(), ngroups, dbarg, maskarg, src2arg);
 533         else
 534             k.args(srcarg, src.cols, (int)src.total(), ngroups, dbarg, maskarg);
 535     }
 536     else
 537     {
 538         if (haveSrc2)
 539             k.args(srcarg, src.cols, (int)src.total(), ngroups, dbarg, src2arg);
 540         else
 541             k.args(srcarg, src.cols, (int)src.total(), ngroups, dbarg);
 542     }
 543
 544     size_t globalsize = ngroups * wgs;
 545     if (k.run(1, &globalsize, &wgs, false))
 546     {
 547         typedef Scalar (*part_sum)(Mat m);
 548         part_sum funcs[3] = { ocl_part_sum<int>, ocl_part_sum<float>, ocl_part_sum<double> },
 549                 func = funcs[ddepth - CV_32S];
 550
 551         Mat mres = db.getMat(ACCESS_READ);
 552         if (calc2)
 553             const_cast<Scalar &>(res2) = func(mres.colRange(ngroups, dbsize));
 554
 555         res = func(mres.colRange(0, ngroups));
 556         return true;
 557     }
 558     return false;
 559 }
 560
 561 #endif
 562
 563 }
 564
 565 cv::Scalar cv::sum( InputArray _src )
 566 {
 567 #ifdef HAVE_OPENCL
 568     Scalar _res;
 569     CV_OCL_RUN_(_src.isUMat() && _src.dims() <= 2,
 570                 ocl_sum(_src, _res, OCL_OP_SUM),
 571                 _res)
 572 #endif
 573
 574     Mat src = _src.getMat();
 575     int k, cn = src.channels(), depth = src.depth();
 576
 577 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
 578     size_t total_size = src.total();
 579     int rows = src.size[0], cols = (int)(total_size/rows);
 580     if( src.dims == 2 || (src.isContinuous() && cols > 0 && (size_t)rows*cols == total_size) )
 581     {
 582         IppiSize sz = { cols, rows };
 583         int type = src.type();
 584         typedef IppStatus (CV_STDCALL* ippiSumFuncHint)(const void*, int, IppiSize, double *, IppHintAlgorithm);
 585         typedef IppStatus (CV_STDCALL* ippiSumFuncNoHint)(const void*, int, IppiSize, double *);
 586         ippiSumFuncHint ippFuncHint =
 587             type == CV_32FC1 ? (ippiSumFuncHint)ippiSum_32f_C1R :
 588             type == CV_32FC3 ? (ippiSumFuncHint)ippiSum_32f_C3R :
 589             type == CV_32FC4 ? (ippiSumFuncHint)ippiSum_32f_C4R :
 590             0;
 591         ippiSumFuncNoHint ippFuncNoHint =
 592             type == CV_8UC1 ? (ippiSumFuncNoHint)ippiSum_8u_C1R :
 593             type == CV_8UC3 ? (ippiSumFuncNoHint)ippiSum_8u_C3R :
 594             type == CV_8UC4 ? (ippiSumFuncNoHint)ippiSum_8u_C4R :
 595             type == CV_16UC1 ? (ippiSumFuncNoHint)ippiSum_16u_C1R :
 596             type == CV_16UC3 ? (ippiSumFuncNoHint)ippiSum_16u_C3R :
 597             type == CV_16UC4 ? (ippiSumFuncNoHint)ippiSum_16u_C4R :
 598             type == CV_16SC1 ? (ippiSumFuncNoHint)ippiSum_16s_C1R :
 599             type == CV_16SC3 ? (ippiSumFuncNoHint)ippiSum_16s_C3R :
 600             type == CV_16SC4 ? (ippiSumFuncNoHint)ippiSum_16s_C4R :
 601             0;
 602         CV_Assert(!ippFuncHint || !ippFuncNoHint);
 603         if( ippFuncHint || ippFuncNoHint )
 604         {
 605             Ipp64f res[4];
 606             IppStatus ret = ippFuncHint ? ippFuncHint(src.data, (int)src.step[0], sz, res, ippAlgHintAccurate) :
 607                             ippFuncNoHint(src.data, (int)src.step[0], sz, res);
 608             if( ret >= 0 )
 609             {
 610                 Scalar sc;
 611                 for( int i = 0; i < cn; i++ )
 612                     sc[i] = res[i];
 613                 return sc;
 614             }
 615             setIppErrorStatus();
 616         }
 617     }
 618 #endif
 619
 620     SumFunc func = getSumFunc(depth);
 621
 622     CV_Assert( cn <= 4 && func != 0 );
 623
 624     const Mat* arrays[] = {&src, 0};
 625     uchar* ptrs[1];
 626     NAryMatIterator it(arrays, ptrs);
 627     Scalar s;
 628     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
 629     int j, count = 0;
 630     AutoBuffer<int> _buf;
 631     int* buf = (int*)&s[0];
 632     size_t esz = 0;
 633     bool blockSum = depth < CV_32S;
 634
 635     if( blockSum )
 636     {
 637         intSumBlockSize = depth <= CV_8S ? (1 << 23) : (1 << 15);
 638         blockSize = std::min(blockSize, intSumBlockSize);
 639         _buf.allocate(cn);
 640         buf = _buf;
 641
 642         for( k = 0; k < cn; k++ )
 643             buf[k] = 0;
 644         esz = src.elemSize();
 645     }
 646
 647     for( size_t i = 0; i < it.nplanes; i++, ++it )
 648     {
 649         for( j = 0; j < total; j += blockSize )
 650         {
 651             int bsz = std::min(total - j, blockSize);
 652             func( ptrs[0], 0, (uchar*)buf, bsz, cn );
 653             count += bsz;
 654             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
 655             {
 656                 for( k = 0; k < cn; k++ )
 657                 {
 658                     s[k] += buf[k];
 659                     buf[k] = 0;
 660                 }
 661                 count = 0;
 662             }
 663             ptrs[0] += bsz*esz;
 664         }
 665     }
 666     return s;
 667 }
 668
 669 #ifdef HAVE_OPENCL
 670
 671 namespace cv {
 672
 673 static bool ocl_countNonZero( InputArray _src, int & res )
 674 {
 675     int type = _src.type(), depth = CV_MAT_DEPTH(type), kercn = ocl::predictOptimalVectorWidth(_src);
 676     bool doubleSupport = ocl::Device::getDefault().doubleFPConfig() > 0;
 677
 678     if (depth == CV_64F && !doubleSupport)
 679         return false;
 680
 681     int dbsize = ocl::Device::getDefault().maxComputeUnits();
 682     size_t wgs = ocl::Device::getDefault().maxWorkGroupSize();
 683
 684     int wgs2_aligned = 1;
 685     while (wgs2_aligned < (int)wgs)
 686         wgs2_aligned <<= 1;
 687     wgs2_aligned >>= 1;
 688
 689     ocl::Kernel k("reduce", ocl::core::reduce_oclsrc,
 690                   format("-D srcT=%s -D srcT1=%s -D cn=1 -D OP_COUNT_NON_ZERO"
 691                          " -D WGS=%d -D kercn=%d -D WGS2_ALIGNED=%d%s%s",
 692                          ocl::typeToStr(CV_MAKE_TYPE(depth, kercn)),
 693                          ocl::typeToStr(depth), (int)wgs, kercn,
 694                          wgs2_aligned, doubleSupport ? " -D DOUBLE_SUPPORT" : "",
 695                          _src.isContinuous() ? " -D HAVE_SRC_CONT" : ""));
 696     if (k.empty())
 697         return false;
 698
 699     UMat src = _src.getUMat(), db(1, dbsize, CV_32SC1);
 700     k.args(ocl::KernelArg::ReadOnlyNoSize(src), src.cols, (int)src.total(),
 701            dbsize, ocl::KernelArg::PtrWriteOnly(db));
 702
 703     size_t globalsize = dbsize * wgs;
 704     if (k.run(1, &globalsize, &wgs, true))
 705         return res = saturate_cast<int>(cv::sum(db.getMat(ACCESS_READ))[0]), true;
 706     return false;
 707 }
 708
 709 }
 710
 711 #endif
 712
 713 int cv::countNonZero( InputArray _src )
 714 {
 715     int type = _src.type(), cn = CV_MAT_CN(type);
 716     CV_Assert( cn == 1 );
 717
 718 #ifdef HAVE_OPENCL
 719     int res = -1;
 720     CV_OCL_RUN_(_src.isUMat() && _src.dims() <= 2,
 721                 ocl_countNonZero(_src, res),
 722                 res)
 723 #endif
 724
 725     Mat src = _src.getMat();
 726
 727 #if defined HAVE_IPP && !defined HAVE_IPP_ICV_ONLY && 0
 728     if (src.dims <= 2 || src.isContinuous())
 729     {
 730         IppiSize roiSize = { src.cols, src.rows };
 731         Ipp32s count = 0, srcstep = (Ipp32s)src.step;
 732         IppStatus status = (IppStatus)-1;
 733
 734         if (src.isContinuous())
 735         {
 736             roiSize.width = (Ipp32s)src.total();
 737             roiSize.height = 1;
 738             srcstep = (Ipp32s)src.total() * CV_ELEM_SIZE(type);
 739         }
 740
 741         int depth = CV_MAT_DEPTH(type);
 742         if (depth == CV_8U)
 743             status = ippiCountInRange_8u_C1R((const Ipp8u *)src.data, srcstep, roiSize, &count, 0, 0);
 744         else if (depth == CV_32F)
 745             status = ippiCountInRange_32f_C1R((const Ipp32f *)src.data, srcstep, roiSize, &count, 0, 0);
 746
 747         if (status >= 0)
 748             return (Ipp32s)src.total() - count;
 749         setIppErrorStatus();
 750     }
 751 #endif
 752
 753     CountNonZeroFunc func = getCountNonZeroTab(src.depth());
 754     CV_Assert( func != 0 );
 755
 756     const Mat* arrays[] = {&src, 0};
 757     uchar* ptrs[1];
 758     NAryMatIterator it(arrays, ptrs);
 759     int total = (int)it.size, nz = 0;
 760
 761     for( size_t i = 0; i < it.nplanes; i++, ++it )
 762         nz += func( ptrs[0], total );
 763
 764     return nz;
 765 }
 766
 767 cv::Scalar cv::mean( InputArray _src, InputArray _mask )
 768 {
 769     Mat src = _src.getMat(), mask = _mask.getMat();
 770     CV_Assert( mask.empty() || mask.type() == CV_8U );
 771
 772     int k, cn = src.channels(), depth = src.depth();
 773
 774 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
 775     size_t total_size = src.total();
 776     int rows = src.size[0], cols = (int)(total_size/rows);
 777     if( src.dims == 2 || (src.isContinuous() && mask.isContinuous() && cols > 0 && (size_t)rows*cols == total_size) )
 778     {
 779         IppiSize sz = { cols, rows };
 780         int type = src.type();
 781         if( !mask.empty() )
 782         {
 783             typedef IppStatus (CV_STDCALL* ippiMaskMeanFuncC1)(const void *, int, void *, int, IppiSize, Ipp64f *);
 784             ippiMaskMeanFuncC1 ippFuncC1 =
 785             type == CV_8UC1 ? (ippiMaskMeanFuncC1)ippiMean_8u_C1MR :
 786             type == CV_16UC1 ? (ippiMaskMeanFuncC1)ippiMean_16u_C1MR :
 787             type == CV_32FC1 ? (ippiMaskMeanFuncC1)ippiMean_32f_C1MR :
 788             0;
 789             if( ippFuncC1 )
 790             {
 791                 Ipp64f res;
 792                 if( ippFuncC1(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, &res) >= 0 )
 793                     return Scalar(res);
 794                 setIppErrorStatus();
 795             }
 796             typedef IppStatus (CV_STDCALL* ippiMaskMeanFuncC3)(const void *, int, void *, int, IppiSize, int, Ipp64f *);
 797             ippiMaskMeanFuncC3 ippFuncC3 =
 798             type == CV_8UC3 ? (ippiMaskMeanFuncC3)ippiMean_8u_C3CMR :
 799             type == CV_16UC3 ? (ippiMaskMeanFuncC3)ippiMean_16u_C3CMR :
 800             type == CV_32FC3 ? (ippiMaskMeanFuncC3)ippiMean_32f_C3CMR :
 801             0;
 802             if( ippFuncC3 )
 803             {
 804                 Ipp64f res1, res2, res3;
 805                 if( ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 1, &res1) >= 0 &&
 806                     ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 2, &res2) >= 0 &&
 807                     ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 3, &res3) >= 0 )
 808                 {
 809                     return Scalar(res1, res2, res3);
 810                 }
 811                 setIppErrorStatus();
 812             }
 813         }
 814         else
 815         {
 816             typedef IppStatus (CV_STDCALL* ippiMeanFuncHint)(const void*, int, IppiSize, double *, IppHintAlgorithm);
 817             typedef IppStatus (CV_STDCALL* ippiMeanFuncNoHint)(const void*, int, IppiSize, double *);
 818             ippiMeanFuncHint ippFuncHint =
 819                 type == CV_32FC1 ? (ippiMeanFuncHint)ippiMean_32f_C1R :
 820                 type == CV_32FC3 ? (ippiMeanFuncHint)ippiMean_32f_C3R :
 821                 type == CV_32FC4 ? (ippiMeanFuncHint)ippiMean_32f_C4R :
 822                 0;
 823             ippiMeanFuncNoHint ippFuncNoHint =
 824                 type == CV_8UC1 ? (ippiMeanFuncNoHint)ippiMean_8u_C1R :
 825                 type == CV_8UC3 ? (ippiMeanFuncNoHint)ippiMean_8u_C3R :
 826                 type == CV_8UC4 ? (ippiMeanFuncNoHint)ippiMean_8u_C4R :
 827                 type == CV_16UC1 ? (ippiMeanFuncNoHint)ippiMean_16u_C1R :
 828                 type == CV_16UC3 ? (ippiMeanFuncNoHint)ippiMean_16u_C3R :
 829                 type == CV_16UC4 ? (ippiMeanFuncNoHint)ippiMean_16u_C4R :
 830                 type == CV_16SC1 ? (ippiMeanFuncNoHint)ippiMean_16s_C1R :
 831                 type == CV_16SC3 ? (ippiMeanFuncNoHint)ippiMean_16s_C3R :
 832                 type == CV_16SC4 ? (ippiMeanFuncNoHint)ippiMean_16s_C4R :
 833                 0;
 834             // Make sure only zero or one version of the function pointer is valid
 835             CV_Assert(!ippFuncHint || !ippFuncNoHint);
 836             if( ippFuncHint || ippFuncNoHint )
 837             {
 838                 Ipp64f res[4];
 839                 IppStatus ret = ippFuncHint ? ippFuncHint(src.data, (int)src.step[0], sz, res, ippAlgHintAccurate) :
 840                                 ippFuncNoHint(src.data, (int)src.step[0], sz, res);
 841                 if( ret >= 0 )
 842                 {
 843                     Scalar sc;
 844                     for( int i = 0; i < cn; i++ )
 845                         sc[i] = res[i];
 846                     return sc;
 847                 }
 848                 setIppErrorStatus();
 849             }
 850         }
 851     }
 852 #endif
 853
 854     SumFunc func = getSumFunc(depth);
 855
 856     CV_Assert( cn <= 4 && func != 0 );
 857
 858     const Mat* arrays[] = {&src, &mask, 0};
 859     uchar* ptrs[2];
 860     NAryMatIterator it(arrays, ptrs);
 861     Scalar s;
 862     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
 863     int j, count = 0;
 864     AutoBuffer<int> _buf;
 865     int* buf = (int*)&s[0];
 866     bool blockSum = depth <= CV_16S;
 867     size_t esz = 0, nz0 = 0;
 868
 869     if( blockSum )
 870     {
 871         intSumBlockSize = depth <= CV_8S ? (1 << 23) : (1 << 15);
 872         blockSize = std::min(blockSize, intSumBlockSize);
 873         _buf.allocate(cn);
 874         buf = _buf;
 875
 876         for( k = 0; k < cn; k++ )
 877             buf[k] = 0;
 878         esz = src.elemSize();
 879     }
 880
 881     for( size_t i = 0; i < it.nplanes; i++, ++it )
 882     {
 883         for( j = 0; j < total; j += blockSize )
 884         {
 885             int bsz = std::min(total - j, blockSize);
 886             int nz = func( ptrs[0], ptrs[1], (uchar*)buf, bsz, cn );
 887             count += nz;
 888             nz0 += nz;
 889             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
 890             {
 891                 for( k = 0; k < cn; k++ )
 892                 {
 893                     s[k] += buf[k];
 894                     buf[k] = 0;
 895                 }
 896                 count = 0;
 897             }
 898             ptrs[0] += bsz*esz;
 899             if( ptrs[1] )
 900                 ptrs[1] += bsz;
 901         }
 902     }
 903     return s*(nz0 ? 1./nz0 : 0);
 904 }
 905
 906 #ifdef HAVE_OPENCL
 907
 908 namespace cv {
 909
 910 static bool ocl_meanStdDev( InputArray _src, OutputArray _mean, OutputArray _sdv, InputArray _mask )
 911 {
 912     bool haveMask = _mask.kind() != _InputArray::NONE;
 913     int nz = haveMask ? -1 : (int)_src.total();
 914     Scalar mean, stddev;
 915
 916     {
 917         int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
 918         bool doubleSupport = ocl::Device::getDefault().doubleFPConfig() > 0,
 919                 isContinuous = _src.isContinuous();
 920         int groups = ocl::Device::getDefault().maxComputeUnits();
 921         size_t wgs = ocl::Device::getDefault().maxWorkGroupSize();
 922
 923         int ddepth = std::max(CV_32S, depth), sqddepth = std::max(CV_32F, depth),
 924                 dtype = CV_MAKE_TYPE(ddepth, cn),
 925                 sqdtype = CV_MAKETYPE(sqddepth, cn);
 926         CV_Assert(!haveMask || _mask.type() == CV_8UC1);
 927
 928         int wgs2_aligned = 1;
 929         while (wgs2_aligned < (int)wgs)
 930             wgs2_aligned <<= 1;
 931         wgs2_aligned >>= 1;
 932
 933         if ( (!doubleSupport && depth == CV_64F) || cn > 4 )
 934             return false;
 935
 936         char cvt[2][40];
 937         String opts = format("-D srcT=%s -D srcT1=%s -D dstT=%s -D dstT1=%s -D sqddepth=%d"
 938                              " -D sqdstT=%s -D sqdstT1=%s -D convertToSDT=%s -D cn=%d%s"
 939                              " -D convertToDT=%s -D WGS=%d -D WGS2_ALIGNED=%d%s%s",
 940                              ocl::typeToStr(type), ocl::typeToStr(depth),
 941                              ocl::typeToStr(dtype), ocl::typeToStr(ddepth), sqddepth,
 942                              ocl::typeToStr(sqdtype), ocl::typeToStr(sqddepth),
 943                              ocl::convertTypeStr(depth, sqddepth, cn, cvt[0]),
 944                              cn, isContinuous ? " -D HAVE_SRC_CONT" : "",
 945                              ocl::convertTypeStr(depth, ddepth, cn, cvt[1]),
 946                              (int)wgs, wgs2_aligned, haveMask ? " -D HAVE_MASK" : "",
 947                              doubleSupport ? " -D DOUBLE_SUPPORT" : "");
 948
 949         ocl::Kernel k("meanStdDev", ocl::core::meanstddev_oclsrc, opts);
 950         if (k.empty())
 951             return false;
 952
 953         int dbsize = groups * ((haveMask ? CV_ELEM_SIZE1(CV_32S) : 0) +
 954                                CV_ELEM_SIZE(sqdtype) + CV_ELEM_SIZE(dtype));
 955         UMat src = _src.getUMat(), db(1, dbsize, CV_8UC1), mask = _mask.getUMat();
 956
 957         ocl::KernelArg srcarg = ocl::KernelArg::ReadOnlyNoSize(src),
 958                 dbarg = ocl::KernelArg::PtrWriteOnly(db),
 959                 maskarg = ocl::KernelArg::ReadOnlyNoSize(mask);
 960
 961         if (haveMask)
 962             k.args(srcarg, src.cols, (int)src.total(), groups, dbarg, maskarg);
 963         else
 964             k.args(srcarg, src.cols, (int)src.total(), groups, dbarg);
 965
 966         size_t globalsize = groups * wgs;
 967         if (!k.run(1, &globalsize, &wgs, false))
 968             return false;
 969
 970         typedef Scalar (* part_sum)(Mat m);
 971         part_sum funcs[3] = { ocl_part_sum<int>, ocl_part_sum<float>, ocl_part_sum<double> };
 972         Mat dbm = db.getMat(ACCESS_READ);
 973
 974         mean = funcs[ddepth - CV_32S](Mat(1, groups, dtype, dbm.data));
 975         stddev = funcs[sqddepth - CV_32S](Mat(1, groups, sqdtype, dbm.data + groups * CV_ELEM_SIZE(dtype)));
 976
 977         if (haveMask)
 978             nz = saturate_cast<int>(funcs[0](Mat(1, groups, CV_32SC1, dbm.data +
 979                                                  groups * (CV_ELEM_SIZE(dtype) +
 980                                                            CV_ELEM_SIZE(sqdtype))))[0]);
 981     }
 982
 983     double total = nz != 0 ? 1.0 / nz : 0;
 984     int k, j, cn = _src.channels();
 985     for (int i = 0; i < cn; ++i)
 986     {
 987         mean[i] *= total;
 988         stddev[i] = std::sqrt(std::max(stddev[i] * total - mean[i] * mean[i] , 0.));
 989     }
 990
 991     for( j = 0; j < 2; j++ )
 992     {
 993         const double * const sptr = j == 0 ? &mean[0] : &stddev[0];
 994         _OutputArray _dst = j == 0 ? _mean : _sdv;
 995         if( !_dst.needed() )
 996             continue;
 997
 998         if( !_dst.fixedSize() )
 999             _dst.create(cn, 1, CV_64F, -1, true);
1000         Mat dst = _dst.getMat();
1001         int dcn = (int)dst.total();
1002         CV_Assert( dst.type() == CV_64F && dst.isContinuous() &&
1003                    (dst.cols == 1 || dst.rows == 1) && dcn >= cn );
1004         double* dptr = dst.ptr<double>();
1005         for( k = 0; k < cn; k++ )
1006             dptr[k] = sptr[k];
1007         for( ; k < dcn; k++ )
1008             dptr[k] = 0;
1009     }
1010
1011     return true;
1012 }
1013
1014 }
1015
1016 #endif
1017
1018 void cv::meanStdDev( InputArray _src, OutputArray _mean, OutputArray _sdv, InputArray _mask )
1019 {
1020     CV_OCL_RUN(_src.isUMat() && _src.dims() <= 2,
1021                ocl_meanStdDev(_src, _mean, _sdv, _mask))
1022
1023     Mat src = _src.getMat(), mask = _mask.getMat();
1024     CV_Assert( mask.empty() || mask.type() == CV_8UC1 );
1025
1026     int k, cn = src.channels(), depth = src.depth();
1027
1028 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
1029     size_t total_size = src.total();
1030     int rows = src.size[0], cols = (int)(total_size/rows);
1031     if( src.dims == 2 || (src.isContinuous() && mask.isContinuous() && cols > 0 && (size_t)rows*cols == total_size) )
1032     {
1033         Ipp64f mean_temp[3];
1034         Ipp64f stddev_temp[3];
1035         Ipp64f *pmean = &mean_temp[0];
1036         Ipp64f *pstddev = &stddev_temp[0];
1037         Mat mean, stddev;
1038         int dcn_mean = -1;
1039         if( _mean.needed() )
1040         {
1041             if( !_mean.fixedSize() )
1042                 _mean.create(cn, 1, CV_64F, -1, true);
1043             mean = _mean.getMat();
1044             dcn_mean = (int)mean.total();
1045             pmean = (Ipp64f *)mean.data;
1046         }
1047         int dcn_stddev = -1;
1048         if( _sdv.needed() )
1049         {
1050             if( !_sdv.fixedSize() )
1051                 _sdv.create(cn, 1, CV_64F, -1, true);
1052             stddev = _sdv.getMat();
1053             dcn_stddev = (int)stddev.total();
1054             pstddev = (Ipp64f *)stddev.data;
1055         }
1056         for( int c = cn; c < dcn_mean; c++ )
1057             pmean[c] = 0;
1058         for( int c = cn; c < dcn_stddev; c++ )
1059             pstddev[c] = 0;
1060         IppiSize sz = { cols, rows };
1061         int type = src.type();
1062         if( !mask.empty() )
1063         {
1064             typedef IppStatus (CV_STDCALL* ippiMaskMeanStdDevFuncC1)(const void *, int, void *, int, IppiSize, Ipp64f *, Ipp64f *);
1065             ippiMaskMeanStdDevFuncC1 ippFuncC1 =
1066             type == CV_8UC1 ? (ippiMaskMeanStdDevFuncC1)ippiMean_StdDev_8u_C1MR :
1067             type == CV_16UC1 ? (ippiMaskMeanStdDevFuncC1)ippiMean_StdDev_16u_C1MR :
1068             type == CV_32FC1 ? (ippiMaskMeanStdDevFuncC1)ippiMean_StdDev_32f_C1MR :
1069             0;
1070             if( ippFuncC1 )
1071             {
1072                 if( ippFuncC1(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, pmean, pstddev) >= 0 )
1073                     return;
1074                 setIppErrorStatus();
1075             }
1076             typedef IppStatus (CV_STDCALL* ippiMaskMeanStdDevFuncC3)(const void *, int, void *, int, IppiSize, int, Ipp64f *, Ipp64f *);
1077             ippiMaskMeanStdDevFuncC3 ippFuncC3 =
1078             type == CV_8UC3 ? (ippiMaskMeanStdDevFuncC3)ippiMean_StdDev_8u_C3CMR :
1079             type == CV_16UC3 ? (ippiMaskMeanStdDevFuncC3)ippiMean_StdDev_16u_C3CMR :
1080             type == CV_32FC3 ? (ippiMaskMeanStdDevFuncC3)ippiMean_StdDev_32f_C3CMR :
1081             0;
1082             if( ippFuncC3 )
1083             {
1084                 if( ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 1, &pmean[0], &pstddev[0]) >= 0 &&
1085                     ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 2, &pmean[1], &pstddev[1]) >= 0 &&
1086                     ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 3, &pmean[2], &pstddev[2]) >= 0 )
1087                     return;
1088                 setIppErrorStatus();
1089             }
1090         }
1091         else
1092         {
1093             typedef IppStatus (CV_STDCALL* ippiMeanStdDevFuncC1)(const void *, int, IppiSize, Ipp64f *, Ipp64f *);
1094             ippiMeanStdDevFuncC1 ippFuncC1 =
1095             type == CV_8UC1 ? (ippiMeanStdDevFuncC1)ippiMean_StdDev_8u_C1R :
1096             type == CV_16UC1 ? (ippiMeanStdDevFuncC1)ippiMean_StdDev_16u_C1R :
1097 #if (IPP_VERSION_X100 >= 801)
1098             type == CV_32FC1 ? (ippiMeanStdDevFuncC1)ippiMean_StdDev_32f_C1R ://Aug 2013: bug in IPP 7.1, 8.0
1099 #endif
1100             0;
1101             if( ippFuncC1 )
1102             {
1103                 if( ippFuncC1(src.data, (int)src.step[0], sz, pmean, pstddev) >= 0 )
1104                     return;
1105                 setIppErrorStatus();
1106             }
1107             typedef IppStatus (CV_STDCALL* ippiMeanStdDevFuncC3)(const void *, int, IppiSize, int, Ipp64f *, Ipp64f *);
1108             ippiMeanStdDevFuncC3 ippFuncC3 =
1109             type == CV_8UC3 ? (ippiMeanStdDevFuncC3)ippiMean_StdDev_8u_C3CR :
1110             type == CV_16UC3 ? (ippiMeanStdDevFuncC3)ippiMean_StdDev_16u_C3CR :
1111             type == CV_32FC3 ? (ippiMeanStdDevFuncC3)ippiMean_StdDev_32f_C3CR :
1112             0;
1113             if( ippFuncC3 )
1114             {
1115                 if( ippFuncC3(src.data, (int)src.step[0], sz, 1, &pmean[0], &pstddev[0]) >= 0 &&
1116                     ippFuncC3(src.data, (int)src.step[0], sz, 2, &pmean[1], &pstddev[1]) >= 0 &&
1117                     ippFuncC3(src.data, (int)src.step[0], sz, 3, &pmean[2], &pstddev[2]) >= 0 )
1118                     return;
1119                 setIppErrorStatus();
1120             }
1121         }
1122     }
1123 #endif
1124
1125
1126     SumSqrFunc func = getSumSqrTab(depth);
1127
1128     CV_Assert( func != 0 );
1129
1130     const Mat* arrays[] = {&src, &mask, 0};
1131     uchar* ptrs[2];
1132     NAryMatIterator it(arrays, ptrs);
1133     int total = (int)it.size, blockSize = total, intSumBlockSize = 0;
1134     int j, count = 0, nz0 = 0;
1135     AutoBuffer<double> _buf(cn*4);
1136     double *s = (double*)_buf, *sq = s + cn;
1137     int *sbuf = (int*)s, *sqbuf = (int*)sq;
1138     bool blockSum = depth <= CV_16S, blockSqSum = depth <= CV_8S;
1139     size_t esz = 0;
1140
1141     for( k = 0; k < cn; k++ )
1142         s[k] = sq[k] = 0;
1143
1144     if( blockSum )
1145     {
1146         intSumBlockSize = 1 << 15;
1147         blockSize = std::min(blockSize, intSumBlockSize);
1148         sbuf = (int*)(sq + cn);
1149         if( blockSqSum )
1150             sqbuf = sbuf + cn;
1151         for( k = 0; k < cn; k++ )
1152             sbuf[k] = sqbuf[k] = 0;
1153         esz = src.elemSize();
1154     }
1155
1156     for( size_t i = 0; i < it.nplanes; i++, ++it )
1157     {
1158         for( j = 0; j < total; j += blockSize )
1159         {
1160             int bsz = std::min(total - j, blockSize);
1161             int nz = func( ptrs[0], ptrs[1], (uchar*)sbuf, (uchar*)sqbuf, bsz, cn );
1162             count += nz;
1163             nz0 += nz;
1164             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
1165             {
1166                 for( k = 0; k < cn; k++ )
1167                 {
1168                     s[k] += sbuf[k];
1169                     sbuf[k] = 0;
1170                 }
1171                 if( blockSqSum )
1172                 {
1173                     for( k = 0; k < cn; k++ )
1174                     {
1175                         sq[k] += sqbuf[k];
1176                         sqbuf[k] = 0;
1177                     }
1178                 }
1179                 count = 0;
1180             }
1181             ptrs[0] += bsz*esz;
1182             if( ptrs[1] )
1183                 ptrs[1] += bsz;
1184         }
1185     }
1186
1187     double scale = nz0 ? 1./nz0 : 0.;
1188     for( k = 0; k < cn; k++ )
1189     {
1190         s[k] *= scale;
1191         sq[k] = std::sqrt(std::max(sq[k]*scale - s[k]*s[k], 0.));
1192     }
1193
1194     for( j = 0; j < 2; j++ )
1195     {
1196         const double* sptr = j == 0 ? s : sq;
1197         _OutputArray _dst = j == 0 ? _mean : _sdv;
1198         if( !_dst.needed() )
1199             continue;
1200
1201         if( !_dst.fixedSize() )
1202             _dst.create(cn, 1, CV_64F, -1, true);
1203         Mat dst = _dst.getMat();
1204         int dcn = (int)dst.total();
1205         CV_Assert( dst.type() == CV_64F && dst.isContinuous() &&
1206                    (dst.cols == 1 || dst.rows == 1) && dcn >= cn );
1207         double* dptr = dst.ptr<double>();
1208         for( k = 0; k < cn; k++ )
1209             dptr[k] = sptr[k];
1210         for( ; k < dcn; k++ )
1211             dptr[k] = 0;
1212     }
1213 }
1214
1215 /****************************************************************************************\
1216 *                                       minMaxLoc                                        *
1217 \****************************************************************************************/
1218
1219 namespace cv
1220 {
1221
1222 template<typename T, typename WT> static void
1223 minMaxIdx_( const T* src, const uchar* mask, WT* _minVal, WT* _maxVal,
1224             size_t* _minIdx, size_t* _maxIdx, int len, size_t startIdx )
1225 {
1226     WT minVal = *_minVal, maxVal = *_maxVal;
1227     size_t minIdx = *_minIdx, maxIdx = *_maxIdx;
1228
1229     if( !mask )
1230     {
1231         for( int i = 0; i < len; i++ )
1232         {
1233             T val = src[i];
1234             if( val < minVal )
1235             {
1236                 minVal = val;
1237                 minIdx = startIdx + i;
1238             }
1239             if( val > maxVal )
1240             {
1241                 maxVal = val;
1242                 maxIdx = startIdx + i;
1243             }
1244         }
1245     }
1246     else
1247     {
1248         for( int i = 0; i < len; i++ )
1249         {
1250             T val = src[i];
1251             if( mask[i] && val < minVal )
1252             {
1253                 minVal = val;
1254                 minIdx = startIdx + i;
1255             }
1256             if( mask[i] && val > maxVal )
1257             {
1258                 maxVal = val;
1259                 maxIdx = startIdx + i;
1260             }
1261         }
1262     }
1263
1264     *_minIdx = minIdx;
1265     *_maxIdx = maxIdx;
1266     *_minVal = minVal;
1267     *_maxVal = maxVal;
1268 }
1269
1270 static void minMaxIdx_8u(const uchar* src, const uchar* mask, int* minval, int* maxval,
1271                          size_t* minidx, size_t* maxidx, int len, size_t startidx )
1272 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
1273
1274 static void minMaxIdx_8s(const schar* src, const uchar* mask, int* minval, int* maxval,
1275                          size_t* minidx, size_t* maxidx, int len, size_t startidx )
1276 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
1277
1278 static void minMaxIdx_16u(const ushort* src, const uchar* mask, int* minval, int* maxval,
1279                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
1280 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
1281
1282 static void minMaxIdx_16s(const short* src, const uchar* mask, int* minval, int* maxval,
1283                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
1284 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
1285
1286 static void minMaxIdx_32s(const int* src, const uchar* mask, int* minval, int* maxval,
1287                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
1288 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
1289
1290 static void minMaxIdx_32f(const float* src, const uchar* mask, float* minval, float* maxval,
1291                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
1292 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
1293
1294 static void minMaxIdx_64f(const double* src, const uchar* mask, double* minval, double* maxval,
1295                           size_t* minidx, size_t* maxidx, int len, size_t startidx )
1296 { minMaxIdx_(src, mask, minval, maxval, minidx, maxidx, len, startidx ); }
1297
1298 typedef void (*MinMaxIdxFunc)(const uchar*, const uchar*, int*, int*, size_t*, size_t*, int, size_t);
1299
1300 static MinMaxIdxFunc getMinmaxTab(int depth)
1301 {
1302     static MinMaxIdxFunc minmaxTab[] =
1303     {
1304         (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_8u), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_8s),
1305         (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_16u), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_16s),
1306         (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_32s),
1307         (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_32f), (MinMaxIdxFunc)GET_OPTIMIZED(minMaxIdx_64f),
1308         0
1309     };
1310
1311     return minmaxTab[depth];
1312 }
1313
1314 static void ofs2idx(const Mat& a, size_t ofs, int* idx)
1315 {
1316     int i, d = a.dims;
1317     if( ofs > 0 )
1318     {
1319         ofs--;
1320         for( i = d-1; i >= 0; i-- )
1321         {
1322             int sz = a.size[i];
1323             idx[i] = (int)(ofs % sz);
1324             ofs /= sz;
1325         }
1326     }
1327     else
1328     {
1329         for( i = d-1; i >= 0; i-- )
1330             idx[i] = -1;
1331     }
1332 }
1333
1334 #ifdef HAVE_OPENCL
1335
1336 template <typename T>
1337 void getMinMaxRes(const Mat & db, double * minVal, double * maxVal,
1338                   int* minLoc, int* maxLoc,
1339                   int groupnum, int cols, double * maxVal2)
1340 {
1341     uint index_max = std::numeric_limits<uint>::max();
1342     T minval = std::numeric_limits<T>::max();
1343     T maxval = std::numeric_limits<T>::min() > 0 ? -std::numeric_limits<T>::max() : std::numeric_limits<T>::min(), maxval2 = maxval;
1344     uint minloc = index_max, maxloc = index_max;
1345
1346     int index = 0;
1347     const T * minptr = NULL, * maxptr = NULL, * maxptr2 = NULL;
1348     const uint * minlocptr = NULL, * maxlocptr = NULL;
1349     if (minVal || minLoc)
1350     {
1351         minptr = (const T *)db.data;
1352         index += sizeof(T) * groupnum;
1353     }
1354     if (maxVal || maxLoc)
1355     {
1356         maxptr = (const T *)(db.data + index);
1357         index += sizeof(T) * groupnum;
1358     }
1359     if (minLoc)
1360     {
1361         minlocptr = (uint *)(db.data + index);
1362         index += sizeof(uint) * groupnum;
1363     }
1364     if (maxLoc)
1365     {
1366         maxlocptr = (uint *)(db.data + index);
1367         index += sizeof(uint) * groupnum;
1368     }
1369     if (maxVal2)
1370         maxptr2 = (const T *)(db.data + index);
1371
1372     for (int i = 0; i < groupnum; i++)
1373     {
1374         if (minptr && minptr[i] <= minval)
1375         {
1376             if (minptr[i] == minval)
1377             {
1378                 if (minlocptr)
1379                     minloc = std::min(minlocptr[i], minloc);
1380             }
1381             else
1382             {
1383                 if (minlocptr)
1384                     minloc = minlocptr[i];
1385                 minval = minptr[i];
1386             }
1387         }
1388         if (maxptr && maxptr[i] >= maxval)
1389         {
1390             if (maxptr[i] == maxval)
1391             {
1392                 if (maxlocptr)
1393                     maxloc = std::min(maxlocptr[i], maxloc);
1394             }
1395             else
1396             {
1397                 if (maxlocptr)
1398                     maxloc = maxlocptr[i];
1399                 maxval = maxptr[i];
1400             }
1401         }
1402         if (maxptr2 && maxptr2[i] > maxval2)
1403             maxval2 = maxptr2[i];
1404     }
1405     bool zero_mask = (minLoc && minloc == index_max) ||
1406             (maxLoc && maxloc == index_max);
1407
1408     if (minVal)
1409         *minVal = zero_mask ? 0 : (double)minval;
1410     if (maxVal)
1411         *maxVal = zero_mask ? 0 : (double)maxval;
1412     if (maxVal2)
1413         *maxVal2 = zero_mask ? 0 : (double)maxval2;
1414
1415     if (minLoc)
1416     {
1417         minLoc[0] = zero_mask ? -1 : minloc / cols;
1418         minLoc[1] = zero_mask ? -1 : minloc % cols;
1419     }
1420     if (maxLoc)
1421     {
1422         maxLoc[0] = zero_mask ? -1 : maxloc / cols;
1423         maxLoc[1] = zero_mask ? -1 : maxloc % cols;
1424     }
1425 }
1426
1427 typedef void (*getMinMaxResFunc)(const Mat & db, double * minVal, double * maxVal,
1428                                  int * minLoc, int *maxLoc, int gropunum, int cols, double * maxVal2);
1429
1430 static bool ocl_minMaxIdx( InputArray _src, double* minVal, double* maxVal, int* minLoc, int* maxLoc, InputArray _mask,
1431                            int ddepth = -1, bool absValues = false, InputArray _src2 = noArray(), double * maxVal2 = NULL)
1432 {
1433     const ocl::Device & dev = ocl::Device::getDefault();
1434     bool doubleSupport = dev.doubleFPConfig() > 0, haveMask = !_mask.empty(),
1435         haveSrc2 = _src2.kind() != _InputArray::NONE;
1436     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type),
1437             kercn = haveMask ? cn : std::min(4, ocl::predictOptimalVectorWidth(_src));
1438
1439     CV_Assert( (cn == 1 && (!haveMask || _mask.type() == CV_8U)) ||
1440               (cn >= 1 && (!haveMask || haveSrc2) && !minLoc && !maxLoc) );
1441
1442     if (ddepth < 0)
1443         ddepth = depth;
1444
1445     CV_Assert(!haveSrc2 || _src2.type() == type);
1446
1447     if ((depth == CV_64F || ddepth == CV_64F) && !doubleSupport)
1448         return false;
1449
1450     int groupnum = dev.maxComputeUnits();
1451     size_t wgs = dev.maxWorkGroupSize();
1452
1453     int wgs2_aligned = 1;
1454     while (wgs2_aligned < (int)wgs)
1455         wgs2_aligned <<= 1;
1456     wgs2_aligned >>= 1;
1457
1458     bool needMinVal = minVal || minLoc, needMinLoc = minLoc != NULL,
1459             needMaxVal = maxVal || maxLoc, needMaxLoc = maxLoc != NULL;
1460
1461     // in case of mask we must know whether mask is filled with zeros or not
1462     // so let's calculate min or max location, if it's undefined, so mask is zeros
1463     if (!(needMaxLoc || needMinLoc) && haveMask)
1464     {
1465         if (needMinVal)
1466             needMinLoc = true;
1467         else
1468             needMaxVal = true;
1469     }
1470
1471     char cvt[40];
1472     String opts = format("-D DEPTH_%d -D srcT1=%s%s -D WGS=%d -D srcT=%s"
1473                          " -D WGS2_ALIGNED=%d%s%s%s -D kercn=%d%s%s%s%s"
1474                          " -D dstT1=%s -D dstT=%s -D convertToDT=%s%s%s%s%s",
1475                          depth, ocl::typeToStr(depth), haveMask ? " -D HAVE_MASK" : "", (int)wgs,
1476                          ocl::typeToStr(CV_MAKE_TYPE(depth, kercn)), wgs2_aligned,
1477                          doubleSupport ? " -D DOUBLE_SUPPORT" : "",
1478                          _src.isContinuous() ? " -D HAVE_SRC_CONT" : "",
1479                          _mask.isContinuous() ? " -D HAVE_MASK_CONT" : "", kercn,
1480                          needMinVal ? " -D NEED_MINVAL" : "", needMaxVal ? " -D NEED_MAXVAL" : "",
1481                          needMinLoc ? " -D NEED_MINLOC" : "", needMaxLoc ? " -D NEED_MAXLOC" : "",
1482                          ocl::typeToStr(ddepth), ocl::typeToStr(CV_MAKE_TYPE(ddepth, kercn)),
1483                          ocl::convertTypeStr(depth, ddepth, kercn, cvt), absValues ? " -D OP_ABS" : "",
1484                          haveSrc2 ? " -D HAVE_SRC2" : "", maxVal2 ? " -D OP_CALC2" : "",
1485                          haveSrc2 && _src2.isContinuous() ? " -D HAVE_SRC2_CONT" : "");
1486
1487     printf("%s\n", opts.c_str());
1488
1489     ocl::Kernel k("minmaxloc", ocl::core::minmaxloc_oclsrc, opts);
1490     if (k.empty())
1491         return false;
1492
1493     int esz = CV_ELEM_SIZE(ddepth), esz32s = CV_ELEM_SIZE1(CV_32S),
1494             dbsize = groupnum * ((needMinVal ? esz : 0) + (needMaxVal ? esz : 0) +
1495                                  (needMinLoc ? esz32s : 0) + (needMaxLoc ? esz32s : 0) +
1496                                  (maxVal2 ? esz : 0));
1497     UMat src = _src.getUMat(), src2 = _src2.getUMat(), db(1, dbsize, CV_8UC1), mask = _mask.getUMat();
1498
1499     if (cn > 1)
1500     {
1501         src = src.reshape(1);
1502         src2 = src2.reshape(1);
1503     }
1504
1505     if (haveSrc2)
1506     {
1507         if (!haveMask)
1508             k.args(ocl::KernelArg::ReadOnlyNoSize(src), src.cols, (int)src.total(),
1509                    groupnum, ocl::KernelArg::PtrWriteOnly(db), ocl::KernelArg::ReadOnlyNoSize(src2));
1510         else
1511             k.args(ocl::KernelArg::ReadOnlyNoSize(src), src.cols, (int)src.total(),
1512                    groupnum, ocl::KernelArg::PtrWriteOnly(db), ocl::KernelArg::ReadOnlyNoSize(mask),
1513                    ocl::KernelArg::ReadOnlyNoSize(src2));
1514     }
1515     else
1516     {
1517         if (!haveMask)
1518             k.args(ocl::KernelArg::ReadOnlyNoSize(src), src.cols, (int)src.total(),
1519                    groupnum, ocl::KernelArg::PtrWriteOnly(db));
1520         else
1521             k.args(ocl::KernelArg::ReadOnlyNoSize(src), src.cols, (int)src.total(),
1522                    groupnum, ocl::KernelArg::PtrWriteOnly(db), ocl::KernelArg::ReadOnlyNoSize(mask));
1523     }
1524
1525     size_t globalsize = groupnum * wgs;
1526     if (!k.run(1, &globalsize, &wgs, false))
1527         return false;
1528
1529     static const getMinMaxResFunc functab[7] =
1530     {
1531         getMinMaxRes<uchar>,
1532         getMinMaxRes<char>,
1533         getMinMaxRes<ushort>,
1534         getMinMaxRes<short>,
1535         getMinMaxRes<int>,
1536         getMinMaxRes<float>,
1537         getMinMaxRes<double>
1538     };
1539
1540     getMinMaxResFunc func = functab[ddepth];
1541
1542     int locTemp[2];
1543     func(db.getMat(ACCESS_READ), minVal, maxVal,
1544          needMinLoc ? minLoc ? minLoc : locTemp : minLoc,
1545          needMaxLoc ? maxLoc ? maxLoc : locTemp : maxLoc,
1546          groupnum, src.cols, maxVal2);
1547
1548     return true;
1549 }
1550
1551 #endif
1552
1553 }
1554
1555 void cv::minMaxIdx(InputArray _src, double* minVal,
1556                    double* maxVal, int* minIdx, int* maxIdx,
1557                    InputArray _mask)
1558 {
1559     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
1560     CV_Assert( (cn == 1 && (_mask.empty() || _mask.type() == CV_8U)) ||
1561         (cn > 1 && _mask.empty() && !minIdx && !maxIdx) );
1562
1563     CV_OCL_RUN(_src.isUMat() && _src.dims() <= 2  && (_mask.empty() || _src.size() == _mask.size()),
1564                ocl_minMaxIdx(_src, minVal, maxVal, minIdx, maxIdx, _mask))
1565
1566     Mat src = _src.getMat(), mask = _mask.getMat();
1567
1568 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
1569     size_t total_size = src.total();
1570     int rows = src.size[0], cols = (int)(total_size/rows);
1571     if( src.dims == 2 || (src.isContinuous() && mask.isContinuous() && cols > 0 && (size_t)rows*cols == total_size) )
1572     {
1573         IppiSize sz = { cols * cn, rows };
1574
1575         if( !mask.empty() )
1576         {
1577             typedef IppStatus (CV_STDCALL* ippiMaskMinMaxIndxFuncC1)(const void *, int, const void *, int,
1578                                                                      IppiSize, Ipp32f *, Ipp32f *, IppiPoint *, IppiPoint *);
1579
1580             CV_SUPPRESS_DEPRECATED_START
1581             ippiMaskMinMaxIndxFuncC1 ippFuncC1 =
1582                 type == CV_8UC1 ? (ippiMaskMinMaxIndxFuncC1)ippiMinMaxIndx_8u_C1MR :
1583                 type == CV_8SC1 ? (ippiMaskMinMaxIndxFuncC1)ippiMinMaxIndx_8s_C1MR :
1584                 type == CV_16UC1 ? (ippiMaskMinMaxIndxFuncC1)ippiMinMaxIndx_16u_C1MR :
1585                 type == CV_32FC1 ? (ippiMaskMinMaxIndxFuncC1)ippiMinMaxIndx_32f_C1MR : 0;
1586             CV_SUPPRESS_DEPRECATED_END
1587
1588             if( ippFuncC1 )
1589             {
1590                 Ipp32f min, max;
1591                 IppiPoint minp, maxp;
1592                 if( ippFuncC1(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, &min, &max, &minp, &maxp) >= 0 )
1593                 {
1594                     if( minVal )
1595                         *minVal = (double)min;
1596                     if( maxVal )
1597                         *maxVal = (double)max;
1598                     if( !minp.x && !minp.y && !maxp.x && !maxp.y && !mask.data[0] )
1599                         minp.x = maxp.x = -1;
1600                     if( minIdx )
1601                     {
1602                         size_t minidx = minp.y * cols + minp.x + 1;
1603                         ofs2idx(src, minidx, minIdx);
1604                     }
1605                     if( maxIdx )
1606                     {
1607                         size_t maxidx = maxp.y * cols + maxp.x + 1;
1608                         ofs2idx(src, maxidx, maxIdx);
1609                     }
1610                     return;
1611                 }
1612                 setIppErrorStatus();
1613             }
1614         }
1615         else
1616         {
1617             typedef IppStatus (CV_STDCALL* ippiMinMaxIndxFuncC1)(const void *, int, IppiSize, Ipp32f *, Ipp32f *, IppiPoint *, IppiPoint *);
1618
1619             CV_SUPPRESS_DEPRECATED_START
1620             ippiMinMaxIndxFuncC1 ippFuncC1 =
1621                 depth == CV_8U ? (ippiMinMaxIndxFuncC1)ippiMinMaxIndx_8u_C1R :
1622                 depth == CV_8S ? (ippiMinMaxIndxFuncC1)ippiMinMaxIndx_8s_C1R :
1623                 depth == CV_16U ? (ippiMinMaxIndxFuncC1)ippiMinMaxIndx_16u_C1R :
1624                 depth == CV_32F ? (ippiMinMaxIndxFuncC1)ippiMinMaxIndx_32f_C1R : 0;
1625             CV_SUPPRESS_DEPRECATED_END
1626
1627             if( ippFuncC1 )
1628             {
1629                 Ipp32f min, max;
1630                 IppiPoint minp, maxp;
1631                 if( ippFuncC1(src.data, (int)src.step[0], sz, &min, &max, &minp, &maxp) >= 0 )
1632                 {
1633                     if( minVal )
1634                         *minVal = (double)min;
1635                     if( maxVal )
1636                         *maxVal = (double)max;
1637                     if( minIdx )
1638                     {
1639                         size_t minidx = minp.y * cols + minp.x + 1;
1640                         ofs2idx(src, minidx, minIdx);
1641                     }
1642                     if( maxIdx )
1643                     {
1644                         size_t maxidx = maxp.y * cols + maxp.x + 1;
1645                         ofs2idx(src, maxidx, maxIdx);
1646                     }
1647                     return;
1648                 }
1649                 setIppErrorStatus();
1650             }
1651         }
1652     }
1653 #endif
1654
1655     MinMaxIdxFunc func = getMinmaxTab(depth);
1656     CV_Assert( func != 0 );
1657
1658     const Mat* arrays[] = {&src, &mask, 0};
1659     uchar* ptrs[2];
1660     NAryMatIterator it(arrays, ptrs);
1661
1662     size_t minidx = 0, maxidx = 0;
1663     int iminval = INT_MAX, imaxval = INT_MIN;
1664     float fminval = FLT_MAX, fmaxval = -FLT_MAX;
1665     double dminval = DBL_MAX, dmaxval = -DBL_MAX;
1666     size_t startidx = 1;
1667     int *minval = &iminval, *maxval = &imaxval;
1668     int planeSize = (int)it.size*cn;
1669
1670     if( depth == CV_32F )
1671         minval = (int*)&fminval, maxval = (int*)&fmaxval;
1672     else if( depth == CV_64F )
1673         minval = (int*)&dminval, maxval = (int*)&dmaxval;
1674
1675     for( size_t i = 0; i < it.nplanes; i++, ++it, startidx += planeSize )
1676         func( ptrs[0], ptrs[1], minval, maxval, &minidx, &maxidx, planeSize, startidx );
1677
1678     if( minidx == 0 )
1679         dminval = dmaxval = 0;
1680     else if( depth == CV_32F )
1681         dminval = fminval, dmaxval = fmaxval;
1682     else if( depth <= CV_32S )
1683         dminval = iminval, dmaxval = imaxval;
1684
1685     if( minVal )
1686         *minVal = dminval;
1687     if( maxVal )
1688         *maxVal = dmaxval;
1689
1690     if( minIdx )
1691         ofs2idx(src, minidx, minIdx);
1692     if( maxIdx )
1693         ofs2idx(src, maxidx, maxIdx);
1694 }
1695
1696 void cv::minMaxLoc( InputArray _img, double* minVal, double* maxVal,
1697                     Point* minLoc, Point* maxLoc, InputArray mask )
1698 {
1699     CV_Assert(_img.dims() <= 2);
1700
1701     minMaxIdx(_img, minVal, maxVal, (int*)minLoc, (int*)maxLoc, mask);
1702     if( minLoc )
1703         std::swap(minLoc->x, minLoc->y);
1704     if( maxLoc )
1705         std::swap(maxLoc->x, maxLoc->y);
1706 }
1707
1708 /****************************************************************************************\
1709 *                                         norm                                           *
1710 \****************************************************************************************/
1711
1712 namespace cv
1713 {
1714
1715 float normL2Sqr_(const float* a, const float* b, int n)
1716 {
1717     int j = 0; float d = 0.f;
1718 #if CV_SSE
1719     if( USE_SSE2 )
1720     {
1721         float CV_DECL_ALIGNED(16) buf[4];
1722         __m128 d0 = _mm_setzero_ps(), d1 = _mm_setzero_ps();
1723
1724         for( ; j <= n - 8; j += 8 )
1725         {
1726             __m128 t0 = _mm_sub_ps(_mm_loadu_ps(a + j), _mm_loadu_ps(b + j));
1727             __m128 t1 = _mm_sub_ps(_mm_loadu_ps(a + j + 4), _mm_loadu_ps(b + j + 4));
1728             d0 = _mm_add_ps(d0, _mm_mul_ps(t0, t0));
1729             d1 = _mm_add_ps(d1, _mm_mul_ps(t1, t1));
1730         }
1731         _mm_store_ps(buf, _mm_add_ps(d0, d1));
1732         d = buf[0] + buf[1] + buf[2] + buf[3];
1733     }
1734     else
1735 #endif
1736     {
1737         for( ; j <= n - 4; j += 4 )
1738         {
1739             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];
1740             d += t0*t0 + t1*t1 + t2*t2 + t3*t3;
1741         }
1742     }
1743
1744     for( ; j < n; j++ )
1745     {
1746         float t = a[j] - b[j];
1747         d += t*t;
1748     }
1749     return d;
1750 }
1751
1752
1753 float normL1_(const float* a, const float* b, int n)
1754 {
1755     int j = 0; float d = 0.f;
1756 #if CV_SSE
1757     if( USE_SSE2 )
1758     {
1759         float CV_DECL_ALIGNED(16) buf[4];
1760         static const int CV_DECL_ALIGNED(16) absbuf[4] = {0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff};
1761         __m128 d0 = _mm_setzero_ps(), d1 = _mm_setzero_ps();
1762         __m128 absmask = _mm_load_ps((const float*)absbuf);
1763
1764         for( ; j <= n - 8; j += 8 )
1765         {
1766             __m128 t0 = _mm_sub_ps(_mm_loadu_ps(a + j), _mm_loadu_ps(b + j));
1767             __m128 t1 = _mm_sub_ps(_mm_loadu_ps(a + j + 4), _mm_loadu_ps(b + j + 4));
1768             d0 = _mm_add_ps(d0, _mm_and_ps(t0, absmask));
1769             d1 = _mm_add_ps(d1, _mm_and_ps(t1, absmask));
1770         }
1771         _mm_store_ps(buf, _mm_add_ps(d0, d1));
1772         d = buf[0] + buf[1] + buf[2] + buf[3];
1773     }
1774     else
1775 #endif
1776     {
1777         for( ; j <= n - 4; j += 4 )
1778         {
1779             d += std::abs(a[j] - b[j]) + std::abs(a[j+1] - b[j+1]) +
1780                     std::abs(a[j+2] - b[j+2]) + std::abs(a[j+3] - b[j+3]);
1781         }
1782     }
1783
1784     for( ; j < n; j++ )
1785         d += std::abs(a[j] - b[j]);
1786     return d;
1787 }
1788
1789 int normL1_(const uchar* a, const uchar* b, int n)
1790 {
1791     int j = 0, d = 0;
1792 #if CV_SSE
1793     if( USE_SSE2 )
1794     {
1795         __m128i d0 = _mm_setzero_si128();
1796
1797         for( ; j <= n - 16; j += 16 )
1798         {
1799             __m128i t0 = _mm_loadu_si128((const __m128i*)(a + j));
1800             __m128i t1 = _mm_loadu_si128((const __m128i*)(b + j));
1801
1802             d0 = _mm_add_epi32(d0, _mm_sad_epu8(t0, t1));
1803         }
1804
1805         for( ; j <= n - 4; j += 4 )
1806         {
1807             __m128i t0 = _mm_cvtsi32_si128(*(const int*)(a + j));
1808             __m128i t1 = _mm_cvtsi32_si128(*(const int*)(b + j));
1809
1810             d0 = _mm_add_epi32(d0, _mm_sad_epu8(t0, t1));
1811         }
1812         d = _mm_cvtsi128_si32(_mm_add_epi32(d0, _mm_unpackhi_epi64(d0, d0)));
1813     }
1814     else
1815 #endif
1816     {
1817         for( ; j <= n - 4; j += 4 )
1818         {
1819             d += std::abs(a[j] - b[j]) + std::abs(a[j+1] - b[j+1]) +
1820                     std::abs(a[j+2] - b[j+2]) + std::abs(a[j+3] - b[j+3]);
1821         }
1822     }
1823     for( ; j < n; j++ )
1824         d += std::abs(a[j] - b[j]);
1825     return d;
1826 }
1827
1828 static const uchar popCountTable[] =
1829 {
1830     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,
1831     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,
1832     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,
1833     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,
1834     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,
1835     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,
1836     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,
1837     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
1838 };
1839
1840 static const uchar popCountTable2[] =
1841 {
1842     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,
1843     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,
1844     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,
1845     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,
1846     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,
1847     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,
1848     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,
1849     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
1850 };
1851
1852 static const uchar popCountTable4[] =
1853 {
1854     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,
1855     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,
1856     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,
1857     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,
1858     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,
1859     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,
1860     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,
1861     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
1862 };
1863
1864 static int normHamming(const uchar* a, int n)
1865 {
1866     int i = 0, result = 0;
1867 #if CV_NEON
1868     {
1869         uint32x4_t bits = vmovq_n_u32(0);
1870         for (; i <= n - 16; i += 16) {
1871             uint8x16_t A_vec = vld1q_u8 (a + i);
1872             uint8x16_t bitsSet = vcntq_u8 (A_vec);
1873             uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
1874             uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
1875             bits = vaddq_u32(bits, bitSet4);
1876         }
1877         uint64x2_t bitSet2 = vpaddlq_u32 (bits);
1878         result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
1879         result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
1880     }
1881 #endif
1882         for( ; i <= n - 4; i += 4 )
1883             result += popCountTable[a[i]] + popCountTable[a[i+1]] +
1884             popCountTable[a[i+2]] + popCountTable[a[i+3]];
1885     for( ; i < n; i++ )
1886         result += popCountTable[a[i]];
1887     return result;
1888 }
1889
1890 int normHamming(const uchar* a, const uchar* b, int n)
1891 {
1892     int i = 0, result = 0;
1893 #if CV_NEON
1894     {
1895         uint32x4_t bits = vmovq_n_u32(0);
1896         for (; i <= n - 16; i += 16) {
1897             uint8x16_t A_vec = vld1q_u8 (a + i);
1898             uint8x16_t B_vec = vld1q_u8 (b + i);
1899             uint8x16_t AxorB = veorq_u8 (A_vec, B_vec);
1900             uint8x16_t bitsSet = vcntq_u8 (AxorB);
1901             uint16x8_t bitSet8 = vpaddlq_u8 (bitsSet);
1902             uint32x4_t bitSet4 = vpaddlq_u16 (bitSet8);
1903             bits = vaddq_u32(bits, bitSet4);
1904         }
1905         uint64x2_t bitSet2 = vpaddlq_u32 (bits);
1906         result = vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),0);
1907         result += vgetq_lane_s32 (vreinterpretq_s32_u64(bitSet2),2);
1908     }
1909 #endif
1910         for( ; i <= n - 4; i += 4 )
1911             result += popCountTable[a[i] ^ b[i]] + popCountTable[a[i+1] ^ b[i+1]] +
1912                     popCountTable[a[i+2] ^ b[i+2]] + popCountTable[a[i+3] ^ b[i+3]];
1913     for( ; i < n; i++ )
1914         result += popCountTable[a[i] ^ b[i]];
1915     return result;
1916 }
1917
1918 static int normHamming(const uchar* a, int n, int cellSize)
1919 {
1920     if( cellSize == 1 )
1921         return normHamming(a, n);
1922     const uchar* tab = 0;
1923     if( cellSize == 2 )
1924         tab = popCountTable2;
1925     else if( cellSize == 4 )
1926         tab = popCountTable4;
1927     else
1928         CV_Error( CV_StsBadSize, "bad cell size (not 1, 2 or 4) in normHamming" );
1929     int i = 0, result = 0;
1930 #if CV_ENABLE_UNROLLED
1931     for( ; i <= n - 4; i += 4 )
1932         result += tab[a[i]] + tab[a[i+1]] + tab[a[i+2]] + tab[a[i+3]];
1933 #endif
1934     for( ; i < n; i++ )
1935         result += tab[a[i]];
1936     return result;
1937 }
1938
1939 int normHamming(const uchar* a, const uchar* b, int n, int cellSize)
1940 {
1941     if( cellSize == 1 )
1942         return normHamming(a, b, n);
1943     const uchar* tab = 0;
1944     if( cellSize == 2 )
1945         tab = popCountTable2;
1946     else if( cellSize == 4 )
1947         tab = popCountTable4;
1948     else
1949         CV_Error( CV_StsBadSize, "bad cell size (not 1, 2 or 4) in normHamming" );
1950     int i = 0, result = 0;
1951     #if CV_ENABLE_UNROLLED
1952     for( ; i <= n - 4; i += 4 )
1953         result += tab[a[i] ^ b[i]] + tab[a[i+1] ^ b[i+1]] +
1954                 tab[a[i+2] ^ b[i+2]] + tab[a[i+3] ^ b[i+3]];
1955     #endif
1956     for( ; i < n; i++ )
1957         result += tab[a[i] ^ b[i]];
1958     return result;
1959 }
1960
1961
1962 template<typename T, typename ST> int
1963 normInf_(const T* src, const uchar* mask, ST* _result, int len, int cn)
1964 {
1965     ST result = *_result;
1966     if( !mask )
1967     {
1968         result = std::max(result, normInf<T, ST>(src, len*cn));
1969     }
1970     else
1971     {
1972         for( int i = 0; i < len; i++, src += cn )
1973             if( mask[i] )
1974             {
1975                 for( int k = 0; k < cn; k++ )
1976                     result = std::max(result, ST(std::abs(src[k])));
1977             }
1978     }
1979     *_result = result;
1980     return 0;
1981 }
1982
1983 template<typename T, typename ST> int
1984 normL1_(const T* src, const uchar* mask, ST* _result, int len, int cn)
1985 {
1986     ST result = *_result;
1987     if( !mask )
1988     {
1989         result += normL1<T, ST>(src, len*cn);
1990     }
1991     else
1992     {
1993         for( int i = 0; i < len; i++, src += cn )
1994             if( mask[i] )
1995             {
1996                 for( int k = 0; k < cn; k++ )
1997                     result += std::abs(src[k]);
1998             }
1999     }
2000     *_result = result;
2001     return 0;
2002 }
2003
2004 template<typename T, typename ST> int
2005 normL2_(const T* src, const uchar* mask, ST* _result, int len, int cn)
2006 {
2007     ST result = *_result;
2008     if( !mask )
2009     {
2010         result += normL2Sqr<T, ST>(src, len*cn);
2011     }
2012     else
2013     {
2014         for( int i = 0; i < len; i++, src += cn )
2015             if( mask[i] )
2016             {
2017                 for( int k = 0; k < cn; k++ )
2018                 {
2019                     T v = src[k];
2020                     result += (ST)v*v;
2021                 }
2022             }
2023     }
2024     *_result = result;
2025     return 0;
2026 }
2027
2028 template<typename T, typename ST> int
2029 normDiffInf_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
2030 {
2031     ST result = *_result;
2032     if( !mask )
2033     {
2034         result = std::max(result, normInf<T, ST>(src1, src2, len*cn));
2035     }
2036     else
2037     {
2038         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
2039             if( mask[i] )
2040             {
2041                 for( int k = 0; k < cn; k++ )
2042                     result = std::max(result, (ST)std::abs(src1[k] - src2[k]));
2043             }
2044     }
2045     *_result = result;
2046     return 0;
2047 }
2048
2049 template<typename T, typename ST> int
2050 normDiffL1_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
2051 {
2052     ST result = *_result;
2053     if( !mask )
2054     {
2055         result += normL1<T, ST>(src1, src2, len*cn);
2056     }
2057     else
2058     {
2059         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
2060             if( mask[i] )
2061             {
2062                 for( int k = 0; k < cn; k++ )
2063                     result += std::abs(src1[k] - src2[k]);
2064             }
2065     }
2066     *_result = result;
2067     return 0;
2068 }
2069
2070 template<typename T, typename ST> int
2071 normDiffL2_(const T* src1, const T* src2, const uchar* mask, ST* _result, int len, int cn)
2072 {
2073     ST result = *_result;
2074     if( !mask )
2075     {
2076         result += normL2Sqr<T, ST>(src1, src2, len*cn);
2077     }
2078     else
2079     {
2080         for( int i = 0; i < len; i++, src1 += cn, src2 += cn )
2081             if( mask[i] )
2082             {
2083                 for( int k = 0; k < cn; k++ )
2084                 {
2085                     ST v = src1[k] - src2[k];
2086                     result += v*v;
2087                 }
2088             }
2089     }
2090     *_result = result;
2091     return 0;
2092 }
2093
2094
2095 #define CV_DEF_NORM_FUNC(L, suffix, type, ntype) \
2096     static int norm##L##_##suffix(const type* src, const uchar* mask, ntype* r, int len, int cn) \
2097 { return norm##L##_(src, mask, r, len, cn); } \
2098     static int normDiff##L##_##suffix(const type* src1, const type* src2, \
2099     const uchar* mask, ntype* r, int len, int cn) \
2100 { return normDiff##L##_(src1, src2, mask, r, (int)len, cn); }
2101
2102 #define CV_DEF_NORM_ALL(suffix, type, inftype, l1type, l2type) \
2103     CV_DEF_NORM_FUNC(Inf, suffix, type, inftype) \
2104     CV_DEF_NORM_FUNC(L1, suffix, type, l1type) \
2105     CV_DEF_NORM_FUNC(L2, suffix, type, l2type)
2106
2107 CV_DEF_NORM_ALL(8u, uchar, int, int, int)
2108 CV_DEF_NORM_ALL(8s, schar, int, int, int)
2109 CV_DEF_NORM_ALL(16u, ushort, int, int, double)
2110 CV_DEF_NORM_ALL(16s, short, int, int, double)
2111 CV_DEF_NORM_ALL(32s, int, int, double, double)
2112 CV_DEF_NORM_ALL(32f, float, float, double, double)
2113 CV_DEF_NORM_ALL(64f, double, double, double, double)
2114
2115
2116 typedef int (*NormFunc)(const uchar*, const uchar*, uchar*, int, int);
2117 typedef int (*NormDiffFunc)(const uchar*, const uchar*, const uchar*, uchar*, int, int);
2118
2119 static NormFunc getNormFunc(int normType, int depth)
2120 {
2121     static NormFunc normTab[3][8] =
2122     {
2123         {
2124             (NormFunc)GET_OPTIMIZED(normInf_8u), (NormFunc)GET_OPTIMIZED(normInf_8s), (NormFunc)GET_OPTIMIZED(normInf_16u), (NormFunc)GET_OPTIMIZED(normInf_16s),
2125             (NormFunc)GET_OPTIMIZED(normInf_32s), (NormFunc)GET_OPTIMIZED(normInf_32f), (NormFunc)normInf_64f, 0
2126         },
2127         {
2128             (NormFunc)GET_OPTIMIZED(normL1_8u), (NormFunc)GET_OPTIMIZED(normL1_8s), (NormFunc)GET_OPTIMIZED(normL1_16u), (NormFunc)GET_OPTIMIZED(normL1_16s),
2129             (NormFunc)GET_OPTIMIZED(normL1_32s), (NormFunc)GET_OPTIMIZED(normL1_32f), (NormFunc)normL1_64f, 0
2130         },
2131         {
2132             (NormFunc)GET_OPTIMIZED(normL2_8u), (NormFunc)GET_OPTIMIZED(normL2_8s), (NormFunc)GET_OPTIMIZED(normL2_16u), (NormFunc)GET_OPTIMIZED(normL2_16s),
2133             (NormFunc)GET_OPTIMIZED(normL2_32s), (NormFunc)GET_OPTIMIZED(normL2_32f), (NormFunc)normL2_64f, 0
2134         }
2135     };
2136
2137     return normTab[normType][depth];
2138 }
2139
2140 static NormDiffFunc getNormDiffFunc(int normType, int depth)
2141 {
2142     static NormDiffFunc normDiffTab[3][8] =
2143     {
2144         {
2145             (NormDiffFunc)GET_OPTIMIZED(normDiffInf_8u), (NormDiffFunc)normDiffInf_8s,
2146             (NormDiffFunc)normDiffInf_16u, (NormDiffFunc)normDiffInf_16s,
2147             (NormDiffFunc)normDiffInf_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffInf_32f),
2148             (NormDiffFunc)normDiffInf_64f, 0
2149         },
2150         {
2151             (NormDiffFunc)GET_OPTIMIZED(normDiffL1_8u), (NormDiffFunc)normDiffL1_8s,
2152             (NormDiffFunc)normDiffL1_16u, (NormDiffFunc)normDiffL1_16s,
2153             (NormDiffFunc)normDiffL1_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffL1_32f),
2154             (NormDiffFunc)normDiffL1_64f, 0
2155         },
2156         {
2157             (NormDiffFunc)GET_OPTIMIZED(normDiffL2_8u), (NormDiffFunc)normDiffL2_8s,
2158             (NormDiffFunc)normDiffL2_16u, (NormDiffFunc)normDiffL2_16s,
2159             (NormDiffFunc)normDiffL2_32s, (NormDiffFunc)GET_OPTIMIZED(normDiffL2_32f),
2160             (NormDiffFunc)normDiffL2_64f, 0
2161         }
2162     };
2163
2164     return normDiffTab[normType][depth];
2165 }
2166
2167 #ifdef HAVE_OPENCL
2168
2169 static bool ocl_norm( InputArray _src, int normType, InputArray _mask, double & result )
2170 {
2171     const ocl::Device & d = ocl::Device::getDefault();
2172     int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
2173     bool doubleSupport = d.doubleFPConfig() > 0,
2174             haveMask = _mask.kind() != _InputArray::NONE;
2175
2176     if ( !(normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 || normType == NORM_L2SQR) ||
2177          (!doubleSupport && depth == CV_64F))
2178         return false;
2179
2180     UMat src = _src.getUMat();
2181
2182     if (normType == NORM_INF)
2183     {
2184         if (cn == 1 || !haveMask)
2185             ocl_minMaxIdx(_src, NULL, &result, NULL, NULL, _mask,
2186                           std::max(depth, CV_32S), depth != CV_8U && depth != CV_16U);
2187         else
2188         {
2189             int dbsize = d.maxComputeUnits();
2190             size_t wgs = d.maxWorkGroupSize();
2191
2192             int wgs2_aligned = 1;
2193             while (wgs2_aligned < (int)wgs)
2194                 wgs2_aligned <<= 1;
2195             wgs2_aligned >>= 1;
2196
2197             ocl::Kernel k("reduce", ocl::core::reduce_oclsrc,
2198                           format("-D OP_NORM_INF_MASK -D HAVE_MASK -D DEPTH_%d"
2199                                  " -D srcT=%s -D srcT1=%s -D WGS=%d -D cn=%d -D WGS2_ALIGNED=%d%s%s%s",
2200                                  depth, ocl::typeToStr(type), ocl::typeToStr(depth),
2201                                  wgs, cn, wgs2_aligned, doubleSupport ? " -D DOUBLE_SUPPORT" : "",
2202                                  src.isContinuous() ? " -D HAVE_CONT_SRC" : "",
2203                                  _mask.isContinuous() ? " -D HAVE_MASK_CONT" : ""));
2204             if (k.empty())
2205                 return false;
2206
2207             UMat db(1, dbsize, type), mask = _mask.getUMat();
2208             k.args(ocl::KernelArg::ReadOnlyNoSize(src), src.cols, (int)src.total(),
2209                    dbsize, ocl::KernelArg::PtrWriteOnly(db), ocl::KernelArg::ReadOnlyNoSize(mask));
2210
2211             size_t globalsize = dbsize * wgs;
2212             if (!k.run(1, &globalsize, &wgs, true))
2213                 return false;
2214
2215             minMaxIdx(db.getMat(ACCESS_READ), NULL, &result, NULL, NULL, noArray());
2216         }
2217     }
2218     else if (normType == NORM_L1 || normType == NORM_L2 || normType == NORM_L2SQR)
2219     {
2220         Scalar sc;
2221         bool unstype = depth == CV_8U || depth == CV_16U;
2222
2223         if ( !ocl_sum(haveMask ? src : src.reshape(1), sc, normType == NORM_L2 || normType == NORM_L2SQR ?
2224                     OCL_OP_SUM_SQR : (unstype ? OCL_OP_SUM : OCL_OP_SUM_ABS), _mask) )
2225             return false;
2226
2227         if (!haveMask)
2228             cn = 1;
2229
2230         double s = 0.0;
2231         for (int i = 0; i < cn; ++i)
2232             s += sc[i];
2233
2234         result = normType == NORM_L1 || normType == NORM_L2SQR ? s : std::sqrt(s);
2235     }
2236
2237     return true;
2238 }
2239
2240 #endif
2241
2242 }
2243
2244 double cv::norm( InputArray _src, int normType, InputArray _mask )
2245 {
2246     normType &= NORM_TYPE_MASK;
2247     CV_Assert( normType == NORM_INF || normType == NORM_L1 ||
2248                normType == NORM_L2 || normType == NORM_L2SQR ||
2249                ((normType == NORM_HAMMING || normType == NORM_HAMMING2) && _src.type() == CV_8U) );
2250
2251 #ifdef HAVE_OPENCL
2252     double _result = 0;
2253     CV_OCL_RUN_(_src.isUMat() && _src.dims() <= 2,
2254                 ocl_norm(_src, normType, _mask, _result),
2255                 _result)
2256 #endif
2257
2258     Mat src = _src.getMat(), mask = _mask.getMat();
2259     int depth = src.depth(), cn = src.channels();
2260
2261 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
2262     size_t total_size = src.total();
2263     int rows = src.size[0], cols = (int)(total_size/rows);
2264
2265     if( (src.dims == 2 || (src.isContinuous() && mask.isContinuous()))
2266         && cols > 0 && (size_t)rows*cols == total_size
2267         && (normType == NORM_INF || normType == NORM_L1 ||
2268             normType == NORM_L2 || normType == NORM_L2SQR) )
2269     {
2270         IppiSize sz = { cols, rows };
2271         int type = src.type();
2272         if( !mask.empty() )
2273         {
2274             typedef IppStatus (CV_STDCALL* ippiMaskNormFuncC1)(const void *, int, const void *, int, IppiSize, Ipp64f *);
2275             ippiMaskNormFuncC1 ippFuncC1 =
2276                 normType == NORM_INF ?
2277                 (type == CV_8UC1 ? (ippiMaskNormFuncC1)ippiNorm_Inf_8u_C1MR :
2278                 type == CV_8SC1 ? (ippiMaskNormFuncC1)ippiNorm_Inf_8s_C1MR :
2279 //                type == CV_16UC1 ? (ippiMaskNormFuncC1)ippiNorm_Inf_16u_C1MR :
2280                 type == CV_32FC1 ? (ippiMaskNormFuncC1)ippiNorm_Inf_32f_C1MR :
2281                 0) :
2282             normType == NORM_L1 ?
2283                 (type == CV_8UC1 ? (ippiMaskNormFuncC1)ippiNorm_L1_8u_C1MR :
2284                 type == CV_8SC1 ? (ippiMaskNormFuncC1)ippiNorm_L1_8s_C1MR :
2285                 type == CV_16UC1 ? (ippiMaskNormFuncC1)ippiNorm_L1_16u_C1MR :
2286                 type == CV_32FC1 ? (ippiMaskNormFuncC1)ippiNorm_L1_32f_C1MR :
2287                 0) :
2288             normType == NORM_L2 || normType == NORM_L2SQR ?
2289                 (type == CV_8UC1 ? (ippiMaskNormFuncC1)ippiNorm_L2_8u_C1MR :
2290                 type == CV_8SC1 ? (ippiMaskNormFuncC1)ippiNorm_L2_8s_C1MR :
2291                 type == CV_16UC1 ? (ippiMaskNormFuncC1)ippiNorm_L2_16u_C1MR :
2292                 type == CV_32FC1 ? (ippiMaskNormFuncC1)ippiNorm_L2_32f_C1MR :
2293                 0) : 0;
2294             if( ippFuncC1 )
2295             {
2296                 Ipp64f norm;
2297                 if( ippFuncC1(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, &norm) >= 0 )
2298                     return normType == NORM_L2SQR ? (double)(norm * norm) : (double)norm;
2299
2300                 setIppErrorStatus();
2301             }
2302             typedef IppStatus (CV_STDCALL* ippiMaskNormFuncC3)(const void *, int, const void *, int, IppiSize, int, Ipp64f *);
2303             ippiMaskNormFuncC3 ippFuncC3 =
2304                 normType == NORM_INF ?
2305                 (type == CV_8UC3 ? (ippiMaskNormFuncC3)ippiNorm_Inf_8u_C3CMR :
2306                 type == CV_8SC3 ? (ippiMaskNormFuncC3)ippiNorm_Inf_8s_C3CMR :
2307                 type == CV_16UC3 ? (ippiMaskNormFuncC3)ippiNorm_Inf_16u_C3CMR :
2308                 type == CV_32FC3 ? (ippiMaskNormFuncC3)ippiNorm_Inf_32f_C3CMR :
2309                 0) :
2310             normType == NORM_L1 ?
2311                 (type == CV_8UC3 ? (ippiMaskNormFuncC3)ippiNorm_L1_8u_C3CMR :
2312                 type == CV_8SC3 ? (ippiMaskNormFuncC3)ippiNorm_L1_8s_C3CMR :
2313                 type == CV_16UC3 ? (ippiMaskNormFuncC3)ippiNorm_L1_16u_C3CMR :
2314                 type == CV_32FC3 ? (ippiMaskNormFuncC3)ippiNorm_L1_32f_C3CMR :
2315                 0) :
2316             normType == NORM_L2 || normType == NORM_L2SQR ?
2317                 (type == CV_8UC3 ? (ippiMaskNormFuncC3)ippiNorm_L2_8u_C3CMR :
2318                 type == CV_8SC3 ? (ippiMaskNormFuncC3)ippiNorm_L2_8s_C3CMR :
2319                 type == CV_16UC3 ? (ippiMaskNormFuncC3)ippiNorm_L2_16u_C3CMR :
2320                 type == CV_32FC3 ? (ippiMaskNormFuncC3)ippiNorm_L2_32f_C3CMR :
2321                 0) : 0;
2322             if( ippFuncC3 )
2323             {
2324                 Ipp64f norm1, norm2, norm3;
2325                 if( ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 1, &norm1) >= 0 &&
2326                     ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 2, &norm2) >= 0 &&
2327                     ippFuncC3(src.data, (int)src.step[0], mask.data, (int)mask.step[0], sz, 3, &norm3) >= 0)
2328                 {
2329                     Ipp64f norm =
2330                         normType == NORM_INF ? std::max(std::max(norm1, norm2), norm3) :
2331                         normType == NORM_L1 ? norm1 + norm2 + norm3 :
2332                         normType == NORM_L2 || normType == NORM_L2SQR ? std::sqrt(norm1 * norm1 + norm2 * norm2 + norm3 * norm3) :
2333                         0;
2334                     return normType == NORM_L2SQR ? (double)(norm * norm) : (double)norm;
2335                 }
2336                 setIppErrorStatus();
2337             }
2338         }
2339         else
2340         {
2341             typedef IppStatus (CV_STDCALL* ippiNormFuncHint)(const void *, int, IppiSize, Ipp64f *, IppHintAlgorithm hint);
2342             typedef IppStatus (CV_STDCALL* ippiNormFuncNoHint)(const void *, int, IppiSize, Ipp64f *);
2343             ippiNormFuncHint ippFuncHint =
2344                 normType == NORM_L1 ?
2345                 (type == CV_32FC1 ? (ippiNormFuncHint)ippiNorm_L1_32f_C1R :
2346                 type == CV_32FC3 ? (ippiNormFuncHint)ippiNorm_L1_32f_C3R :
2347                 type == CV_32FC4 ? (ippiNormFuncHint)ippiNorm_L1_32f_C4R :
2348                 0) :
2349                 normType == NORM_L2 || normType == NORM_L2SQR ?
2350                 (type == CV_32FC1 ? (ippiNormFuncHint)ippiNorm_L2_32f_C1R :
2351                 type == CV_32FC3 ? (ippiNormFuncHint)ippiNorm_L2_32f_C3R :
2352                 type == CV_32FC4 ? (ippiNormFuncHint)ippiNorm_L2_32f_C4R :
2353                 0) : 0;
2354             ippiNormFuncNoHint ippFuncNoHint =
2355                 normType == NORM_INF ?
2356                 (type == CV_8UC1 ? (ippiNormFuncNoHint)ippiNorm_Inf_8u_C1R :
2357                 type == CV_8UC3 ? (ippiNormFuncNoHint)ippiNorm_Inf_8u_C3R :
2358                 type == CV_8UC4 ? (ippiNormFuncNoHint)ippiNorm_Inf_8u_C4R :
2359                 type == CV_16UC1 ? (ippiNormFuncNoHint)ippiNorm_Inf_16u_C1R :
2360                 type == CV_16UC3 ? (ippiNormFuncNoHint)ippiNorm_Inf_16u_C3R :
2361                 type == CV_16UC4 ? (ippiNormFuncNoHint)ippiNorm_Inf_16u_C4R :
2362                 type == CV_16SC1 ? (ippiNormFuncNoHint)ippiNorm_Inf_16s_C1R :
2363 #if (IPP_VERSION_X100 >= 801)
2364                 type == CV_16SC3 ? (ippiNormFuncNoHint)ippiNorm_Inf_16s_C3R : //Aug 2013: problem in IPP 7.1, 8.0 : -32768
2365                 type == CV_16SC4 ? (ippiNormFuncNoHint)ippiNorm_Inf_16s_C4R : //Aug 2013: problem in IPP 7.1, 8.0 : -32768
2366 #endif
2367                 type == CV_32FC1 ? (ippiNormFuncNoHint)ippiNorm_Inf_32f_C1R :
2368                 type == CV_32FC3 ? (ippiNormFuncNoHint)ippiNorm_Inf_32f_C3R :
2369                 type == CV_32FC4 ? (ippiNormFuncNoHint)ippiNorm_Inf_32f_C4R :
2370                 0) :
2371                 normType == NORM_L1 ?
2372                 (type == CV_8UC1 ? (ippiNormFuncNoHint)ippiNorm_L1_8u_C1R :
2373                 type == CV_8UC3 ? (ippiNormFuncNoHint)ippiNorm_L1_8u_C3R :
2374                 type == CV_8UC4 ? (ippiNormFuncNoHint)ippiNorm_L1_8u_C4R :
2375                 type == CV_16UC1 ? (ippiNormFuncNoHint)ippiNorm_L1_16u_C1R :
2376                 type == CV_16UC3 ? (ippiNormFuncNoHint)ippiNorm_L1_16u_C3R :
2377                 type == CV_16UC4 ? (ippiNormFuncNoHint)ippiNorm_L1_16u_C4R :
2378                 type == CV_16SC1 ? (ippiNormFuncNoHint)ippiNorm_L1_16s_C1R :
2379                 type == CV_16SC3 ? (ippiNormFuncNoHint)ippiNorm_L1_16s_C3R :
2380                 type == CV_16SC4 ? (ippiNormFuncNoHint)ippiNorm_L1_16s_C4R :
2381                 0) :
2382                 normType == NORM_L2 || normType == NORM_L2SQR ?
2383                 (type == CV_8UC1 ? (ippiNormFuncNoHint)ippiNorm_L2_8u_C1R :
2384                 type == CV_8UC3 ? (ippiNormFuncNoHint)ippiNorm_L2_8u_C3R :
2385                 type == CV_8UC4 ? (ippiNormFuncNoHint)ippiNorm_L2_8u_C4R :
2386                 type == CV_16UC1 ? (ippiNormFuncNoHint)ippiNorm_L2_16u_C1R :
2387                 type == CV_16UC3 ? (ippiNormFuncNoHint)ippiNorm_L2_16u_C3R :
2388                 type == CV_16UC4 ? (ippiNormFuncNoHint)ippiNorm_L2_16u_C4R :
2389                 type == CV_16SC1 ? (ippiNormFuncNoHint)ippiNorm_L2_16s_C1R :
2390                 type == CV_16SC3 ? (ippiNormFuncNoHint)ippiNorm_L2_16s_C3R :
2391                 type == CV_16SC4 ? (ippiNormFuncNoHint)ippiNorm_L2_16s_C4R :
2392                 0) : 0;
2393             // Make sure only zero or one version of the function pointer is valid
2394             CV_Assert(!ippFuncHint || !ippFuncNoHint);
2395             if( ippFuncHint || ippFuncNoHint )
2396             {
2397                 Ipp64f norm_array[4];
2398                 IppStatus ret = ippFuncHint ? ippFuncHint(src.data, (int)src.step[0], sz, norm_array, ippAlgHintAccurate) :
2399                                 ippFuncNoHint(src.data, (int)src.step[0], sz, norm_array);
2400                 if( ret >= 0 )
2401                 {
2402                     Ipp64f norm = (normType == NORM_L2 || normType == NORM_L2SQR) ? norm_array[0] * norm_array[0] : norm_array[0];
2403                     for( int i = 1; i < cn; i++ )
2404                     {
2405                         norm =
2406                             normType == NORM_INF ? std::max(norm, norm_array[i]) :
2407                             normType == NORM_L1 ? norm + norm_array[i] :
2408                             normType == NORM_L2 || normType == NORM_L2SQR ? norm + norm_array[i] * norm_array[i] :
2409                             0;
2410                     }
2411                     return normType == NORM_L2 ? (double)std::sqrt(norm) : (double)norm;
2412                 }
2413                 setIppErrorStatus();
2414             }
2415         }
2416     }
2417 #endif
2418
2419     if( src.isContinuous() && mask.empty() )
2420     {
2421         size_t len = src.total()*cn;
2422         if( len == (size_t)(int)len )
2423         {
2424             if( depth == CV_32F )
2425             {
2426                 const float* data = src.ptr<float>();
2427
2428                 if( normType == NORM_L2 )
2429                 {
2430                     double result = 0;
2431                     GET_OPTIMIZED(normL2_32f)(data, 0, &result, (int)len, 1);
2432                     return std::sqrt(result);
2433                 }
2434                 if( normType == NORM_L2SQR )
2435                 {
2436                     double result = 0;
2437                     GET_OPTIMIZED(normL2_32f)(data, 0, &result, (int)len, 1);
2438                     return result;
2439                 }
2440                 if( normType == NORM_L1 )
2441                 {
2442                     double result = 0;
2443                     GET_OPTIMIZED(normL1_32f)(data, 0, &result, (int)len, 1);
2444                     return result;
2445                 }
2446                 if( normType == NORM_INF )
2447                 {
2448                     float result = 0;
2449                     GET_OPTIMIZED(normInf_32f)(data, 0, &result, (int)len, 1);
2450                     return result;
2451                 }
2452             }
2453             if( depth == CV_8U )
2454             {
2455                 const uchar* data = src.ptr<uchar>();
2456
2457                 if( normType == NORM_HAMMING )
2458                     return normHamming(data, (int)len);
2459
2460                 if( normType == NORM_HAMMING2 )
2461                     return normHamming(data, (int)len, 2);
2462             }
2463         }
2464     }
2465
2466     CV_Assert( mask.empty() || mask.type() == CV_8U );
2467
2468     if( normType == NORM_HAMMING || normType == NORM_HAMMING2 )
2469     {
2470         if( !mask.empty() )
2471         {
2472             Mat temp;
2473             bitwise_and(src, mask, temp);
2474             return norm(temp, normType);
2475         }
2476         int cellSize = normType == NORM_HAMMING ? 1 : 2;
2477
2478         const Mat* arrays[] = {&src, 0};
2479         uchar* ptrs[1];
2480         NAryMatIterator it(arrays, ptrs);
2481         int total = (int)it.size;
2482         int result = 0;
2483
2484         for( size_t i = 0; i < it.nplanes; i++, ++it )
2485             result += normHamming(ptrs[0], total, cellSize);
2486
2487         return result;
2488     }
2489
2490     NormFunc func = getNormFunc(normType >> 1, depth);
2491     CV_Assert( func != 0 );
2492
2493     const Mat* arrays[] = {&src, &mask, 0};
2494     uchar* ptrs[2];
2495     union
2496     {
2497         double d;
2498         int i;
2499         float f;
2500     }
2501     result;
2502     result.d = 0;
2503     NAryMatIterator it(arrays, ptrs);
2504     int j, total = (int)it.size, blockSize = total, intSumBlockSize = 0, count = 0;
2505     bool blockSum = (normType == NORM_L1 && depth <= CV_16S) ||
2506             ((normType == NORM_L2 || normType == NORM_L2SQR) && depth <= CV_8S);
2507     int isum = 0;
2508     int *ibuf = &result.i;
2509     size_t esz = 0;
2510
2511     if( blockSum )
2512     {
2513         intSumBlockSize = (normType == NORM_L1 && depth <= CV_8S ? (1 << 23) : (1 << 15))/cn;
2514         blockSize = std::min(blockSize, intSumBlockSize);
2515         ibuf = &isum;
2516         esz = src.elemSize();
2517     }
2518
2519     for( size_t i = 0; i < it.nplanes; i++, ++it )
2520     {
2521         for( j = 0; j < total; j += blockSize )
2522         {
2523             int bsz = std::min(total - j, blockSize);
2524             func( ptrs[0], ptrs[1], (uchar*)ibuf, bsz, cn );
2525             count += bsz;
2526             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
2527             {
2528                 result.d += isum;
2529                 isum = 0;
2530                 count = 0;
2531             }
2532             ptrs[0] += bsz*esz;
2533             if( ptrs[1] )
2534                 ptrs[1] += bsz;
2535         }
2536     }
2537
2538     if( normType == NORM_INF )
2539     {
2540         if( depth == CV_64F )
2541             ;
2542         else if( depth == CV_32F )
2543             result.d = result.f;
2544         else
2545             result.d = result.i;
2546     }
2547     else if( normType == NORM_L2 )
2548         result.d = std::sqrt(result.d);
2549
2550     return result.d;
2551 }
2552
2553 #ifdef HAVE_OPENCL
2554
2555 namespace cv {
2556
2557 static bool ocl_norm( InputArray _src1, InputArray _src2, int normType, InputArray _mask, double & result )
2558 {
2559     Scalar sc1, sc2;
2560     int type = _src1.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
2561     bool relative = (normType & NORM_RELATIVE) != 0;
2562     normType &= ~NORM_RELATIVE;
2563     bool normsum = normType == NORM_L1 || normType == NORM_L2 || normType == NORM_L2SQR;
2564
2565     if ( !(normType == NORM_INF || normsum) )
2566         return false;
2567
2568     if (normsum)
2569     {
2570         if (!ocl_sum(_src1, sc1, normType == NORM_L2 || normType == NORM_L2SQR ?
2571                      OCL_OP_SUM_SQR : OCL_OP_SUM, _mask, _src2, relative, sc2))
2572             return false;
2573     }
2574     else
2575     {
2576         if (!ocl_minMaxIdx(_src1, NULL, &sc1[0], NULL, NULL, _mask, std::max(CV_32S, depth),
2577                            false, _src2, relative ? &sc2[0] : NULL))
2578             return false;
2579         cn = 1;
2580     }
2581
2582     double s2 = 0;
2583     for (int i = 0; i < cn; ++i)
2584     {
2585         result += sc1[i];
2586         if (relative)
2587             s2 += sc2[i];
2588     }
2589
2590     if (normType == NORM_L2)
2591     {
2592         result = std::sqrt(result);
2593         if (relative)
2594             s2 = std::sqrt(s2);
2595     }
2596
2597     if (relative)
2598         result /= (s2 + DBL_EPSILON);
2599
2600     return true;
2601 }
2602
2603 }
2604
2605 #endif
2606
2607 double cv::norm( InputArray _src1, InputArray _src2, int normType, InputArray _mask )
2608 {
2609     CV_Assert( _src1.sameSize(_src2) && _src1.type() == _src2.type() );
2610
2611 #ifdef HAVE_OPENCL
2612     double _result = 0;
2613     CV_OCL_RUN_(_src1.isUMat(),
2614                 ocl_norm(_src1, _src2, normType, _mask, _result),
2615                 _result)
2616 #endif
2617
2618     if( normType & CV_RELATIVE )
2619     {
2620 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
2621         Mat src1 = _src1.getMat(), src2 = _src2.getMat(), mask = _mask.getMat();
2622
2623         normType &= NORM_TYPE_MASK;
2624         CV_Assert( normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 || normType == NORM_L2SQR ||
2625                 ((normType == NORM_HAMMING || normType == NORM_HAMMING2) && src1.type() == CV_8U) );
2626         size_t total_size = src1.total();
2627         int rows = src1.size[0], cols = (int)(total_size/rows);
2628         if( (src1.dims == 2 || (src1.isContinuous() && src2.isContinuous() && mask.isContinuous()))
2629             && cols > 0 && (size_t)rows*cols == total_size
2630             && (normType == NORM_INF || normType == NORM_L1 ||
2631                 normType == NORM_L2 || normType == NORM_L2SQR) )
2632         {
2633             IppiSize sz = { cols, rows };
2634             int type = src1.type();
2635             if( !mask.empty() )
2636             {
2637                 typedef IppStatus (CV_STDCALL* ippiMaskNormRelFuncC1)(const void *, int, const void *, int, const void *, int, IppiSize, Ipp64f *);
2638                 ippiMaskNormRelFuncC1 ippFuncC1 =
2639                     normType == NORM_INF ?
2640                     (type == CV_8UC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_Inf_8u_C1MR :
2641                     type == CV_8SC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_Inf_8s_C1MR :
2642                     type == CV_16UC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_Inf_16u_C1MR :
2643                     type == CV_32FC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_Inf_32f_C1MR :
2644                     0) :
2645                     normType == NORM_L1 ?
2646                     (type == CV_8UC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_L1_8u_C1MR :
2647                     type == CV_8SC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_L1_8s_C1MR :
2648                     type == CV_16UC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_L1_16u_C1MR :
2649                     type == CV_32FC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_L1_32f_C1MR :
2650                     0) :
2651                     normType == NORM_L2 || normType == NORM_L2SQR ?
2652                     (type == CV_8UC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_L2_8u_C1MR :
2653                     type == CV_8SC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_L2_8s_C1MR :
2654                     type == CV_16UC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_L2_16u_C1MR :
2655                     type == CV_32FC1 ? (ippiMaskNormRelFuncC1)ippiNormRel_L2_32f_C1MR :
2656                     0) : 0;
2657                 if( ippFuncC1 )
2658                 {
2659                     Ipp64f norm;
2660                     if( ippFuncC1(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], mask.data, (int)mask.step[0], sz, &norm) >= 0 )
2661                         return normType == NORM_L2SQR ? (double)(norm * norm) : (double)norm;
2662                     setIppErrorStatus();
2663                 }
2664             }
2665             else
2666             {
2667                 typedef IppStatus (CV_STDCALL* ippiNormRelFuncNoHint)(const void *, int, const void *, int, IppiSize, Ipp64f *);
2668                 typedef IppStatus (CV_STDCALL* ippiNormRelFuncHint)(const void *, int, const void *, int, IppiSize, Ipp64f *, IppHintAlgorithm hint);
2669                 ippiNormRelFuncNoHint ippFuncNoHint =
2670                     normType == NORM_INF ?
2671                     (type == CV_8UC1 ? (ippiNormRelFuncNoHint)ippiNormRel_Inf_8u_C1R :
2672                     type == CV_16UC1 ? (ippiNormRelFuncNoHint)ippiNormRel_Inf_16u_C1R :
2673                     type == CV_16SC1 ? (ippiNormRelFuncNoHint)ippiNormRel_Inf_16s_C1R :
2674                     type == CV_32FC1 ? (ippiNormRelFuncNoHint)ippiNormRel_Inf_32f_C1R :
2675                     0) :
2676                     normType == NORM_L1 ?
2677                     (type == CV_8UC1 ? (ippiNormRelFuncNoHint)ippiNormRel_L1_8u_C1R :
2678                     type == CV_16UC1 ? (ippiNormRelFuncNoHint)ippiNormRel_L1_16u_C1R :
2679                     type == CV_16SC1 ? (ippiNormRelFuncNoHint)ippiNormRel_L1_16s_C1R :
2680                     0) :
2681                     normType == NORM_L2 || normType == NORM_L2SQR ?
2682                     (type == CV_8UC1 ? (ippiNormRelFuncNoHint)ippiNormRel_L2_8u_C1R :
2683                     type == CV_16UC1 ? (ippiNormRelFuncNoHint)ippiNormRel_L2_16u_C1R :
2684                     type == CV_16SC1 ? (ippiNormRelFuncNoHint)ippiNormRel_L2_16s_C1R :
2685                     0) : 0;
2686                 ippiNormRelFuncHint ippFuncHint =
2687                     normType == NORM_L1 ?
2688                     (type == CV_32FC1 ? (ippiNormRelFuncHint)ippiNormRel_L1_32f_C1R :
2689                     0) :
2690                     normType == NORM_L2 || normType == NORM_L2SQR ?
2691                     (type == CV_32FC1 ? (ippiNormRelFuncHint)ippiNormRel_L2_32f_C1R :
2692                     0) : 0;
2693                 if (ippFuncNoHint)
2694                 {
2695                     Ipp64f norm;
2696                     if( ippFuncNoHint(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], sz, &norm) >= 0 )
2697                         return (double)norm;
2698                     setIppErrorStatus();
2699                 }
2700                 if (ippFuncHint)
2701                 {
2702                     Ipp64f norm;
2703                     if( ippFuncHint(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], sz, &norm, ippAlgHintAccurate) >= 0 )
2704                         return (double)norm;
2705                     setIppErrorStatus();
2706                 }
2707             }
2708         }
2709 #endif
2710         return norm(_src1, _src2, normType & ~CV_RELATIVE, _mask)/(norm(_src2, normType, _mask) + DBL_EPSILON);
2711     }
2712
2713     Mat src1 = _src1.getMat(), src2 = _src2.getMat(), mask = _mask.getMat();
2714     int depth = src1.depth(), cn = src1.channels();
2715
2716     normType &= 7;
2717     CV_Assert( normType == NORM_INF || normType == NORM_L1 ||
2718                normType == NORM_L2 || normType == NORM_L2SQR ||
2719               ((normType == NORM_HAMMING || normType == NORM_HAMMING2) && src1.type() == CV_8U) );
2720
2721 #if defined (HAVE_IPP) && (IPP_VERSION_MAJOR >= 7)
2722     size_t total_size = src1.total();
2723     int rows = src1.size[0], cols = (int)(total_size/rows);
2724     if( (src1.dims == 2 || (src1.isContinuous() && src2.isContinuous() && mask.isContinuous()))
2725         && cols > 0 && (size_t)rows*cols == total_size
2726         && (normType == NORM_INF || normType == NORM_L1 ||
2727             normType == NORM_L2 || normType == NORM_L2SQR) )
2728     {
2729         IppiSize sz = { cols, rows };
2730         int type = src1.type();
2731         if( !mask.empty() )
2732         {
2733             typedef IppStatus (CV_STDCALL* ippiMaskNormDiffFuncC1)(const void *, int, const void *, int, const void *, int, IppiSize, Ipp64f *);
2734             ippiMaskNormDiffFuncC1 ippFuncC1 =
2735                 normType == NORM_INF ?
2736                 (type == CV_8UC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_Inf_8u_C1MR :
2737                 type == CV_8SC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_Inf_8s_C1MR :
2738                 type == CV_16UC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_Inf_16u_C1MR :
2739                 type == CV_32FC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_Inf_32f_C1MR :
2740                 0) :
2741                 normType == NORM_L1 ?
2742                 (type == CV_8UC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_L1_8u_C1MR :
2743                 type == CV_8SC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_L1_8s_C1MR :
2744                 type == CV_16UC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_L1_16u_C1MR :
2745                 type == CV_32FC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_L1_32f_C1MR :
2746                 0) :
2747                 normType == NORM_L2 || normType == NORM_L2SQR ?
2748                 (type == CV_8UC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_L2_8u_C1MR :
2749                 type == CV_8SC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_L2_8s_C1MR :
2750                 type == CV_16UC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_L2_16u_C1MR :
2751                 type == CV_32FC1 ? (ippiMaskNormDiffFuncC1)ippiNormDiff_L2_32f_C1MR :
2752                 0) : 0;
2753             if( ippFuncC1 )
2754             {
2755                 Ipp64f norm;
2756                 if( ippFuncC1(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], mask.data, (int)mask.step[0], sz, &norm) >= 0 )
2757                     return normType == NORM_L2SQR ? (double)(norm * norm) : (double)norm;
2758                 setIppErrorStatus();
2759             }
2760             typedef IppStatus (CV_STDCALL* ippiMaskNormDiffFuncC3)(const void *, int, const void *, int, const void *, int, IppiSize, int, Ipp64f *);
2761             ippiMaskNormDiffFuncC3 ippFuncC3 =
2762                 normType == NORM_INF ?
2763                 (type == CV_8UC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_Inf_8u_C3CMR :
2764                 type == CV_8SC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_Inf_8s_C3CMR :
2765                 type == CV_16UC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_Inf_16u_C3CMR :
2766                 type == CV_32FC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_Inf_32f_C3CMR :
2767                 0) :
2768                 normType == NORM_L1 ?
2769                 (type == CV_8UC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_L1_8u_C3CMR :
2770                 type == CV_8SC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_L1_8s_C3CMR :
2771                 type == CV_16UC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_L1_16u_C3CMR :
2772                 type == CV_32FC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_L1_32f_C3CMR :
2773                 0) :
2774                 normType == NORM_L2 || normType == NORM_L2SQR ?
2775                 (type == CV_8UC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_L2_8u_C3CMR :
2776                 type == CV_8SC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_L2_8s_C3CMR :
2777                 type == CV_16UC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_L2_16u_C3CMR :
2778                 type == CV_32FC3 ? (ippiMaskNormDiffFuncC3)ippiNormDiff_L2_32f_C3CMR :
2779                 0) : 0;
2780             if( ippFuncC3 )
2781             {
2782                 Ipp64f norm1, norm2, norm3;
2783                 if( ippFuncC3(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], mask.data, (int)mask.step[0], sz, 1, &norm1) >= 0 &&
2784                     ippFuncC3(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], mask.data, (int)mask.step[0], sz, 2, &norm2) >= 0 &&
2785                     ippFuncC3(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], mask.data, (int)mask.step[0], sz, 3, &norm3) >= 0)
2786                 {
2787                     Ipp64f norm =
2788                         normType == NORM_INF ? std::max(std::max(norm1, norm2), norm3) :
2789                         normType == NORM_L1 ? norm1 + norm2 + norm3 :
2790                         normType == NORM_L2 || normType == NORM_L2SQR ? std::sqrt(norm1 * norm1 + norm2 * norm2 + norm3 * norm3) :
2791                         0;
2792                     return normType == NORM_L2SQR ? (double)(norm * norm) : (double)norm;
2793                 }
2794                 setIppErrorStatus();
2795             }
2796         }
2797         else
2798         {
2799             typedef IppStatus (CV_STDCALL* ippiNormDiffFuncHint)(const void *, int, const void *, int, IppiSize, Ipp64f *, IppHintAlgorithm hint);
2800             typedef IppStatus (CV_STDCALL* ippiNormDiffFuncNoHint)(const void *, int, const void *, int, IppiSize, Ipp64f *);
2801             ippiNormDiffFuncHint ippFuncHint =
2802                 normType == NORM_L1 ?
2803                 (type == CV_32FC1 ? (ippiNormDiffFuncHint)ippiNormDiff_L1_32f_C1R :
2804                 type == CV_32FC3 ? (ippiNormDiffFuncHint)ippiNormDiff_L1_32f_C3R :
2805                 type == CV_32FC4 ? (ippiNormDiffFuncHint)ippiNormDiff_L1_32f_C4R :
2806                 0) :
2807                 normType == NORM_L2 || normType == NORM_L2SQR ?
2808                 (type == CV_32FC1 ? (ippiNormDiffFuncHint)ippiNormDiff_L2_32f_C1R :
2809                 type == CV_32FC3 ? (ippiNormDiffFuncHint)ippiNormDiff_L2_32f_C3R :
2810                 type == CV_32FC4 ? (ippiNormDiffFuncHint)ippiNormDiff_L2_32f_C4R :
2811                 0) : 0;
2812             ippiNormDiffFuncNoHint ippFuncNoHint =
2813                 normType == NORM_INF ?
2814                 (type == CV_8UC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_8u_C1R :
2815                 type == CV_8UC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_8u_C3R :
2816                 type == CV_8UC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_8u_C4R :
2817                 type == CV_16UC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_16u_C1R :
2818                 type == CV_16UC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_16u_C3R :
2819                 type == CV_16UC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_16u_C4R :
2820                 type == CV_16SC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_16s_C1R :
2821 #if (IPP_VERSION_X100 >= 801)
2822                 type == CV_16SC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_16s_C3R : //Aug 2013: problem in IPP 7.1, 8.0 : -32768
2823                 type == CV_16SC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_16s_C4R : //Aug 2013: problem in IPP 7.1, 8.0 : -32768
2824 #endif
2825                 type == CV_32FC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_32f_C1R :
2826                 type == CV_32FC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_32f_C3R :
2827                 type == CV_32FC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_Inf_32f_C4R :
2828                 0) :
2829                 normType == NORM_L1 ?
2830                 (type == CV_8UC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_8u_C1R :
2831                 type == CV_8UC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_8u_C3R :
2832                 type == CV_8UC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_8u_C4R :
2833                 type == CV_16UC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_16u_C1R :
2834                 type == CV_16UC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_16u_C3R :
2835                 type == CV_16UC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_16u_C4R :
2836                 type == CV_16SC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_16s_C1R :
2837                 type == CV_16SC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_16s_C3R :
2838                 type == CV_16SC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L1_16s_C4R :
2839                 0) :
2840                 normType == NORM_L2 || normType == NORM_L2SQR ?
2841                 (type == CV_8UC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_8u_C1R :
2842                 type == CV_8UC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_8u_C3R :
2843                 type == CV_8UC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_8u_C4R :
2844                 type == CV_16UC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_16u_C1R :
2845                 type == CV_16UC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_16u_C3R :
2846                 type == CV_16UC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_16u_C4R :
2847                 type == CV_16SC1 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_16s_C1R :
2848                 type == CV_16SC3 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_16s_C3R :
2849                 type == CV_16SC4 ? (ippiNormDiffFuncNoHint)ippiNormDiff_L2_16s_C4R :
2850                 0) : 0;
2851             // Make sure only zero or one version of the function pointer is valid
2852             CV_Assert(!ippFuncHint || !ippFuncNoHint);
2853             if( ippFuncHint || ippFuncNoHint )
2854             {
2855                 Ipp64f norm_array[4];
2856                 IppStatus ret = ippFuncHint ? ippFuncHint(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], sz, norm_array, ippAlgHintAccurate) :
2857                                 ippFuncNoHint(src1.data, (int)src1.step[0], src2.data, (int)src2.step[0], sz, norm_array);
2858                 if( ret >= 0 )
2859                 {
2860                     Ipp64f norm = (normType == NORM_L2 || normType == NORM_L2SQR) ? norm_array[0] * norm_array[0] : norm_array[0];
2861                     for( int i = 1; i < src1.channels(); i++ )
2862                     {
2863                         norm =
2864                             normType == NORM_INF ? std::max(norm, norm_array[i]) :
2865                             normType == NORM_L1 ? norm + norm_array[i] :
2866                             normType == NORM_L2 || normType == NORM_L2SQR ? norm + norm_array[i] * norm_array[i] :
2867                             0;
2868                     }
2869                     return normType == NORM_L2 ? (double)std::sqrt(norm) : (double)norm;
2870                 }
2871                 setIppErrorStatus();
2872             }
2873         }
2874     }
2875 #endif
2876
2877     if( src1.isContinuous() && src2.isContinuous() && mask.empty() )
2878     {
2879         size_t len = src1.total()*src1.channels();
2880         if( len == (size_t)(int)len )
2881         {
2882             if( src1.depth() == CV_32F )
2883             {
2884                 const float* data1 = src1.ptr<float>();
2885                 const float* data2 = src2.ptr<float>();
2886
2887                 if( normType == NORM_L2 )
2888                 {
2889                     double result = 0;
2890                     GET_OPTIMIZED(normDiffL2_32f)(data1, data2, 0, &result, (int)len, 1);
2891                     return std::sqrt(result);
2892                 }
2893                 if( normType == NORM_L2SQR )
2894                 {
2895                     double result = 0;
2896                     GET_OPTIMIZED(normDiffL2_32f)(data1, data2, 0, &result, (int)len, 1);
2897                     return result;
2898                 }
2899                 if( normType == NORM_L1 )
2900                 {
2901                     double result = 0;
2902                     GET_OPTIMIZED(normDiffL1_32f)(data1, data2, 0, &result, (int)len, 1);
2903                     return result;
2904                 }
2905                 if( normType == NORM_INF )
2906                 {
2907                     float result = 0;
2908                     GET_OPTIMIZED(normDiffInf_32f)(data1, data2, 0, &result, (int)len, 1);
2909                     return result;
2910                 }
2911             }
2912         }
2913     }
2914
2915     CV_Assert( mask.empty() || mask.type() == CV_8U );
2916
2917     if( normType == NORM_HAMMING || normType == NORM_HAMMING2 )
2918     {
2919         if( !mask.empty() )
2920         {
2921             Mat temp;
2922             bitwise_xor(src1, src2, temp);
2923             bitwise_and(temp, mask, temp);
2924             return norm(temp, normType);
2925         }
2926         int cellSize = normType == NORM_HAMMING ? 1 : 2;
2927
2928         const Mat* arrays[] = {&src1, &src2, 0};
2929         uchar* ptrs[2];
2930         NAryMatIterator it(arrays, ptrs);
2931         int total = (int)it.size;
2932         int result = 0;
2933
2934         for( size_t i = 0; i < it.nplanes; i++, ++it )
2935             result += normHamming(ptrs[0], ptrs[1], total, cellSize);
2936
2937         return result;
2938     }
2939
2940     NormDiffFunc func = getNormDiffFunc(normType >> 1, depth);
2941     CV_Assert( func != 0 );
2942
2943     const Mat* arrays[] = {&src1, &src2, &mask, 0};
2944     uchar* ptrs[3];
2945     union
2946     {
2947         double d;
2948         float f;
2949         int i;
2950         unsigned u;
2951     }
2952     result;
2953     result.d = 0;
2954     NAryMatIterator it(arrays, ptrs);
2955     int j, total = (int)it.size, blockSize = total, intSumBlockSize = 0, count = 0;
2956     bool blockSum = (normType == NORM_L1 && depth <= CV_16S) ||
2957             ((normType == NORM_L2 || normType == NORM_L2SQR) && depth <= CV_8S);
2958     unsigned isum = 0;
2959     unsigned *ibuf = &result.u;
2960     size_t esz = 0;
2961
2962     if( blockSum )
2963     {
2964         intSumBlockSize = normType == NORM_L1 && depth <= CV_8S ? (1 << 23) : (1 << 15);
2965         blockSize = std::min(blockSize, intSumBlockSize);
2966         ibuf = &isum;
2967         esz = src1.elemSize();
2968     }
2969
2970     for( size_t i = 0; i < it.nplanes; i++, ++it )
2971     {
2972         for( j = 0; j < total; j += blockSize )
2973         {
2974             int bsz = std::min(total - j, blockSize);
2975             func( ptrs[0], ptrs[1], ptrs[2], (uchar*)ibuf, bsz, cn );
2976             count += bsz;
2977             if( blockSum && (count + blockSize >= intSumBlockSize || (i+1 >= it.nplanes && j+bsz >= total)) )
2978             {
2979                 result.d += isum;
2980                 isum = 0;
2981                 count = 0;
2982             }
2983             ptrs[0] += bsz*esz;
2984             ptrs[1] += bsz*esz;
2985             if( ptrs[2] )
2986                 ptrs[2] += bsz;
2987         }
2988     }
2989
2990     if( normType == NORM_INF )
2991     {
2992         if( depth == CV_64F )
2993             ;
2994         else if( depth == CV_32F )
2995             result.d = result.f;
2996         else
2997             result.d = result.u;
2998     }
2999     else if( normType == NORM_L2 )
3000         result.d = std::sqrt(result.d);
3001
3002     return result.d;
3003 }
3004
3005
3006 ///////////////////////////////////// batch distance ///////////////////////////////////////
3007
3008 namespace cv
3009 {
3010
3011 template<typename _Tp, typename _Rt>
3012 void batchDistL1_(const _Tp* src1, const _Tp* src2, size_t step2,
3013                   int nvecs, int len, _Rt* dist, const uchar* mask)
3014 {
3015     step2 /= sizeof(src2[0]);
3016     if( !mask )
3017     {
3018         for( int i = 0; i < nvecs; i++ )
3019             dist[i] = normL1<_Tp, _Rt>(src1, src2 + step2*i, len);
3020     }
3021     else
3022     {
3023         _Rt val0 = std::numeric_limits<_Rt>::max();
3024         for( int i = 0; i < nvecs; i++ )
3025             dist[i] = mask[i] ? normL1<_Tp, _Rt>(src1, src2 + step2*i, len) : val0;
3026     }
3027 }
3028
3029 template<typename _Tp, typename _Rt>
3030 void batchDistL2Sqr_(const _Tp* src1, const _Tp* src2, size_t step2,
3031                      int nvecs, int len, _Rt* dist, const uchar* mask)
3032 {
3033     step2 /= sizeof(src2[0]);
3034     if( !mask )
3035     {
3036         for( int i = 0; i < nvecs; i++ )
3037             dist[i] = normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len);
3038     }
3039     else
3040     {
3041         _Rt val0 = std::numeric_limits<_Rt>::max();
3042         for( int i = 0; i < nvecs; i++ )
3043             dist[i] = mask[i] ? normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len) : val0;
3044     }
3045 }
3046
3047 template<typename _Tp, typename _Rt>
3048 void batchDistL2_(const _Tp* src1, const _Tp* src2, size_t step2,
3049                   int nvecs, int len, _Rt* dist, const uchar* mask)
3050 {
3051     step2 /= sizeof(src2[0]);
3052     if( !mask )
3053     {
3054         for( int i = 0; i < nvecs; i++ )
3055             dist[i] = std::sqrt(normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len));
3056     }
3057     else
3058     {
3059         _Rt val0 = std::numeric_limits<_Rt>::max();
3060         for( int i = 0; i < nvecs; i++ )
3061             dist[i] = mask[i] ? std::sqrt(normL2Sqr<_Tp, _Rt>(src1, src2 + step2*i, len)) : val0;
3062     }
3063 }
3064
3065 static void batchDistHamming(const uchar* src1, const uchar* src2, size_t step2,
3066                              int nvecs, int len, int* dist, const uchar* mask)
3067 {
3068     step2 /= sizeof(src2[0]);
3069     if( !mask )
3070     {
3071         for( int i = 0; i < nvecs; i++ )
3072             dist[i] = normHamming(src1, src2 + step2*i, len);
3073     }
3074     else
3075     {
3076         int val0 = INT_MAX;
3077         for( int i = 0; i < nvecs; i++ )
3078             dist[i] = mask[i] ? normHamming(src1, src2 + step2*i, len) : val0;
3079     }
3080 }
3081
3082 static void batchDistHamming2(const uchar* src1, const uchar* src2, size_t step2,
3083                               int nvecs, int len, int* dist, const uchar* mask)
3084 {
3085     step2 /= sizeof(src2[0]);
3086     if( !mask )
3087     {
3088         for( int i = 0; i < nvecs; i++ )
3089             dist[i] = normHamming(src1, src2 + step2*i, len, 2);
3090     }
3091     else
3092     {
3093         int val0 = INT_MAX;
3094         for( int i = 0; i < nvecs; i++ )
3095             dist[i] = mask[i] ? normHamming(src1, src2 + step2*i, len, 2) : val0;
3096     }
3097 }
3098
3099 static void batchDistL1_8u32s(const uchar* src1, const uchar* src2, size_t step2,
3100                                int nvecs, int len, int* dist, const uchar* mask)
3101 {
3102     batchDistL1_<uchar, int>(src1, src2, step2, nvecs, len, dist, mask);
3103 }
3104
3105 static void batchDistL1_8u32f(const uchar* src1, const uchar* src2, size_t step2,
3106                                int nvecs, int len, float* dist, const uchar* mask)
3107 {
3108     batchDistL1_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
3109 }
3110
3111 static void batchDistL2Sqr_8u32s(const uchar* src1, const uchar* src2, size_t step2,
3112                                   int nvecs, int len, int* dist, const uchar* mask)
3113 {
3114     batchDistL2Sqr_<uchar, int>(src1, src2, step2, nvecs, len, dist, mask);
3115 }
3116
3117 static void batchDistL2Sqr_8u32f(const uchar* src1, const uchar* src2, size_t step2,
3118                                   int nvecs, int len, float* dist, const uchar* mask)
3119 {
3120     batchDistL2Sqr_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
3121 }
3122
3123 static void batchDistL2_8u32f(const uchar* src1, const uchar* src2, size_t step2,
3124                                int nvecs, int len, float* dist, const uchar* mask)
3125 {
3126     batchDistL2_<uchar, float>(src1, src2, step2, nvecs, len, dist, mask);
3127 }
3128
3129 static void batchDistL1_32f(const float* src1, const float* src2, size_t step2,
3130                              int nvecs, int len, float* dist, const uchar* mask)
3131 {
3132     batchDistL1_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
3133 }
3134
3135 static void batchDistL2Sqr_32f(const float* src1, const float* src2, size_t step2,
3136                                 int nvecs, int len, float* dist, const uchar* mask)
3137 {
3138     batchDistL2Sqr_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
3139 }
3140
3141 static void batchDistL2_32f(const float* src1, const float* src2, size_t step2,
3142                              int nvecs, int len, float* dist, const uchar* mask)
3143 {
3144     batchDistL2_<float, float>(src1, src2, step2, nvecs, len, dist, mask);
3145 }
3146
3147 typedef void (*BatchDistFunc)(const uchar* src1, const uchar* src2, size_t step2,
3148                               int nvecs, int len, uchar* dist, const uchar* mask);
3149
3150
3151 struct BatchDistInvoker : public ParallelLoopBody
3152 {
3153     BatchDistInvoker( const Mat& _src1, const Mat& _src2,
3154                       Mat& _dist, Mat& _nidx, int _K,
3155                       const Mat& _mask, int _update,
3156                       BatchDistFunc _func)
3157     {
3158         src1 = &_src1;
3159         src2 = &_src2;
3160         dist = &_dist;
3161         nidx = &_nidx;
3162         K = _K;
3163         mask = &_mask;
3164         update = _update;
3165         func = _func;
3166     }
3167
3168     void operator()(const Range& range) const
3169     {
3170         AutoBuffer<int> buf(src2->rows);
3171         int* bufptr = buf;
3172
3173         for( int i = range.start; i < range.end; i++ )
3174         {
3175             func(src1->ptr(i), src2->ptr(), src2->step, src2->rows, src2->cols,
3176                  K > 0 ? (uchar*)bufptr : dist->ptr(i), mask->data ? mask->ptr(i) : 0);
3177
3178             if( K > 0 )
3179             {
3180                 int* nidxptr = nidx->ptr<int>(i);
3181                 // since positive float's can be compared just like int's,
3182                 // we handle both CV_32S and CV_32F cases with a single branch
3183                 int* distptr = (int*)dist->ptr(i);
3184
3185                 int j, k;
3186
3187                 for( j = 0; j < src2->rows; j++ )
3188                 {
3189                     int d = bufptr[j];
3190                     if( d < distptr[K-1] )
3191                     {
3192                         for( k = K-2; k >= 0 && distptr[k] > d; k-- )
3193                         {
3194                             nidxptr[k+1] = nidxptr[k];
3195                             distptr[k+1] = distptr[k];
3196                         }
3197                         nidxptr[k+1] = j + update;
3198                         distptr[k+1] = d;
3199                     }
3200                 }
3201             }
3202         }
3203     }
3204
3205     const Mat *src1;
3206     const Mat *src2;
3207     Mat *dist;
3208     Mat *nidx;
3209     const Mat *mask;
3210     int K;
3211     int update;
3212     BatchDistFunc func;
3213 };
3214
3215 }
3216
3217 void cv::batchDistance( InputArray _src1, InputArray _src2,
3218                         OutputArray _dist, int dtype, OutputArray _nidx,
3219                         int normType, int K, InputArray _mask,
3220                         int update, bool crosscheck )
3221 {
3222     Mat src1 = _src1.getMat(), src2 = _src2.getMat(), mask = _mask.getMat();
3223     int type = src1.type();
3224     CV_Assert( type == src2.type() && src1.cols == src2.cols &&
3225                (type == CV_32F || type == CV_8U));
3226     CV_Assert( _nidx.needed() == (K > 0) );
3227
3228     if( dtype == -1 )
3229     {
3230         dtype = normType == NORM_HAMMING || normType == NORM_HAMMING2 ? CV_32S : CV_32F;
3231     }
3232     CV_Assert( (type == CV_8U && dtype == CV_32S) || dtype == CV_32F);
3233
3234     K = std::min(K, src2.rows);
3235
3236     _dist.create(src1.rows, (K > 0 ? K : src2.rows), dtype);
3237     Mat dist = _dist.getMat(), nidx;
3238     if( _nidx.needed() )
3239     {
3240         _nidx.create(dist.size(), CV_32S);
3241         nidx = _nidx.getMat();
3242     }
3243
3244     if( update == 0 && K > 0 )
3245     {
3246         dist = Scalar::all(dtype == CV_32S ? (double)INT_MAX : (double)FLT_MAX);
3247         nidx = Scalar::all(-1);
3248     }
3249
3250     if( crosscheck )
3251     {
3252         CV_Assert( K == 1 && update == 0 && mask.empty() );
3253         Mat tdist, tidx;
3254         batchDistance(src2, src1, tdist, dtype, tidx, normType, K, mask, 0, false);
3255
3256         // if an idx-th element from src1 appeared to be the nearest to i-th element of src2,
3257         // we update the minimum mutual distance between idx-th element of src1 and the whole src2 set.
3258         // As a result, if nidx[idx] = i*, it means that idx-th element of src1 is the nearest
3259         // to i*-th element of src2 and i*-th element of src2 is the closest to idx-th element of src1.
3260         // If nidx[idx] = -1, it means that there is no such ideal couple for it in src2.
3261         // This O(N) procedure is called cross-check and it helps to eliminate some false matches.
3262         if( dtype == CV_32S )
3263         {
3264             for( int i = 0; i < tdist.rows; i++ )
3265             {
3266                 int idx = tidx.at<int>(i);
3267                 int d = tdist.at<int>(i), d0 = dist.at<int>(idx);
3268                 if( d < d0 )
3269                 {
3270                     dist.at<int>(idx) = d;
3271                     nidx.at<int>(idx) = i + update;
3272                 }
3273             }
3274         }
3275         else
3276         {
3277             for( int i = 0; i < tdist.rows; i++ )
3278             {
3279                 int idx = tidx.at<int>(i);
3280                 float d = tdist.at<float>(i), d0 = dist.at<float>(idx);
3281                 if( d < d0 )
3282                 {
3283                     dist.at<float>(idx) = d;
3284                     nidx.at<int>(idx) = i + update;
3285                 }
3286             }
3287         }
3288         return;
3289     }
3290
3291     BatchDistFunc func = 0;
3292     if( type == CV_8U )
3293     {
3294         if( normType == NORM_L1 && dtype == CV_32S )
3295             func = (BatchDistFunc)batchDistL1_8u32s;
3296         else if( normType == NORM_L1 && dtype == CV_32F )
3297             func = (BatchDistFunc)batchDistL1_8u32f;
3298         else if( normType == NORM_L2SQR && dtype == CV_32S )
3299             func = (BatchDistFunc)batchDistL2Sqr_8u32s;
3300         else if( normType == NORM_L2SQR && dtype == CV_32F )
3301             func = (BatchDistFunc)batchDistL2Sqr_8u32f;
3302         else if( normType == NORM_L2 && dtype == CV_32F )
3303             func = (BatchDistFunc)batchDistL2_8u32f;
3304         else if( normType == NORM_HAMMING && dtype == CV_32S )
3305             func = (BatchDistFunc)batchDistHamming;
3306         else if( normType == NORM_HAMMING2 && dtype == CV_32S )
3307             func = (BatchDistFunc)batchDistHamming2;
3308     }
3309     else if( type == CV_32F && dtype == CV_32F )
3310     {
3311         if( normType == NORM_L1 )
3312             func = (BatchDistFunc)batchDistL1_32f;
3313         else if( normType == NORM_L2SQR )
3314             func = (BatchDistFunc)batchDistL2Sqr_32f;
3315         else if( normType == NORM_L2 )
3316             func = (BatchDistFunc)batchDistL2_32f;
3317     }
3318
3319     if( func == 0 )
3320         CV_Error_(CV_StsUnsupportedFormat,
3321                   ("The combination of type=%d, dtype=%d and normType=%d is not supported",
3322                    type, dtype, normType));
3323
3324     parallel_for_(Range(0, src1.rows),
3325                   BatchDistInvoker(src1, src2, dist, nidx, K, mask, update, func));
3326 }
3327
3328
3329 void cv::findNonZero( InputArray _src, OutputArray _idx )
3330 {
3331     Mat src = _src.getMat();
3332     CV_Assert( src.type() == CV_8UC1 );
3333     int n = countNonZero(src);
3334     if( _idx.kind() == _InputArray::MAT && !_idx.getMatRef().isContinuous() )
3335         _idx.release();
3336     _idx.create(n, 1, CV_32SC2);
3337     Mat idx = _idx.getMat();
3338     CV_Assert(idx.isContinuous());
3339     Point* idx_ptr = (Point*)idx.data;
3340
3341     for( int i = 0; i < src.rows; i++ )
3342     {
3343         const uchar* bin_ptr = src.ptr(i);
3344         for( int j = 0; j < src.cols; j++ )
3345             if( bin_ptr[j] )
3346                 *idx_ptr++ = Point(j, i);
3347     }
3348 }
3349
3350 double cv::PSNR(InputArray _src1, InputArray _src2)
3351 {
3352     CV_Assert( _src1.depth() == CV_8U );
3353     double diff = std::sqrt(norm(_src1, _src2, NORM_L2SQR)/(_src1.total()*_src1.channels()));
3354     return 20*log10(255./(diff+DBL_EPSILON));
3355 }
3356
3357
3358 CV_IMPL CvScalar cvSum( const CvArr* srcarr )
3359 {
3360     cv::Scalar sum = cv::sum(cv::cvarrToMat(srcarr, false, true, 1));
3361     if( CV_IS_IMAGE(srcarr) )
3362     {
3363         int coi = cvGetImageCOI((IplImage*)srcarr);
3364         if( coi )
3365         {
3366             CV_Assert( 0 < coi && coi <= 4 );
3367             sum = cv::Scalar(sum[coi-1]);
3368         }
3369     }
3370     return sum;
3371 }
3372
3373 CV_IMPL int cvCountNonZero( const CvArr* imgarr )
3374 {
3375     cv::Mat img = cv::cvarrToMat(imgarr, false, true, 1);
3376     if( img.channels() > 1 )
3377         cv::extractImageCOI(imgarr, img);
3378     return countNonZero(img);
3379 }
3380
3381
3382 CV_IMPL  CvScalar
3383 cvAvg( const void* imgarr, const void* maskarr )
3384 {
3385     cv::Mat img = cv::cvarrToMat(imgarr, false, true, 1);
3386     cv::Scalar mean = !maskarr ? cv::mean(img) : cv::mean(img, cv::cvarrToMat(maskarr));
3387     if( CV_IS_IMAGE(imgarr) )
3388     {
3389         int coi = cvGetImageCOI((IplImage*)imgarr);
3390         if( coi )
3391         {
3392             CV_Assert( 0 < coi && coi <= 4 );
3393             mean = cv::Scalar(mean[coi-1]);
3394         }
3395     }
3396     return mean;
3397 }
3398
3399
3400 CV_IMPL  void
3401 cvAvgSdv( const CvArr* imgarr, CvScalar* _mean, CvScalar* _sdv, const void* maskarr )
3402 {
3403     cv::Scalar mean, sdv;
3404
3405     cv::Mat mask;
3406     if( maskarr )
3407         mask = cv::cvarrToMat(maskarr);
3408
3409     cv::meanStdDev(cv::cvarrToMat(imgarr, false, true, 1), mean, sdv, mask );
3410
3411     if( CV_IS_IMAGE(imgarr) )
3412     {
3413         int coi = cvGetImageCOI((IplImage*)imgarr);
3414         if( coi )
3415         {
3416             CV_Assert( 0 < coi && coi <= 4 );
3417             mean = cv::Scalar(mean[coi-1]);
3418             sdv = cv::Scalar(sdv[coi-1]);
3419         }
3420     }
3421
3422     if( _mean )
3423         *(cv::Scalar*)_mean = mean;
3424     if( _sdv )
3425         *(cv::Scalar*)_sdv = sdv;
3426 }
3427
3428
3429 CV_IMPL void
3430 cvMinMaxLoc( const void* imgarr, double* _minVal, double* _maxVal,
3431              CvPoint* _minLoc, CvPoint* _maxLoc, const void* maskarr )
3432 {
3433     cv::Mat mask, img = cv::cvarrToMat(imgarr, false, true, 1);
3434     if( maskarr )
3435         mask = cv::cvarrToMat(maskarr);
3436     if( img.channels() > 1 )
3437         cv::extractImageCOI(imgarr, img);
3438
3439     cv::minMaxLoc( img, _minVal, _maxVal,
3440                    (cv::Point*)_minLoc, (cv::Point*)_maxLoc, mask );
3441 }
3442
3443
3444 CV_IMPL  double
3445 cvNorm( const void* imgA, const void* imgB, int normType, const void* maskarr )
3446 {
3447     cv::Mat a, mask;
3448     if( !imgA )
3449     {
3450         imgA = imgB;
3451         imgB = 0;
3452     }
3453
3454     a = cv::cvarrToMat(imgA, false, true, 1);
3455     if( maskarr )
3456         mask = cv::cvarrToMat(maskarr);
3457
3458     if( a.channels() > 1 && CV_IS_IMAGE(imgA) && cvGetImageCOI((const IplImage*)imgA) > 0 )
3459         cv::extractImageCOI(imgA, a);
3460
3461     if( !imgB )
3462         return !maskarr ? cv::norm(a, normType) : cv::norm(a, normType, mask);
3463
3464     cv::Mat b = cv::cvarrToMat(imgB, false, true, 1);
3465     if( b.channels() > 1 && CV_IS_IMAGE(imgB) && cvGetImageCOI((const IplImage*)imgB) > 0 )
3466         cv::extractImageCOI(imgB, b);
3467
3468     return !maskarr ? cv::norm(a, b, normType) : cv::norm(a, b, normType, mask);
3469 }