modules/ts/src/ts_func.cpp

   1 #include "precomp.hpp"
   2 #include <float.h>
   3 #include <limits.h>
   4 #include "opencv2/imgproc/types_c.h"
   5
   6 using namespace cv;
   7
   8 namespace cvtest
   9 {
  10
  11 const char* getTypeName( int type )
  12 {
  13     static const char* type_names[] = { "8u", "8s", "16u", "16s", "32s", "32f", "64f", "ptr" };
  14     return type_names[CV_MAT_DEPTH(type)];
  15 }
  16
  17 int typeByName( const char* name )
  18 {
  19     int i;
  20     for( i = 0; i < CV_DEPTH_MAX; i++ )
  21         if( strcmp(name, getTypeName(i)) == 0 )
  22             return i;
  23     return -1;
  24 }
  25
  26 string vec2str( const string& sep, const int* v, size_t nelems )
  27 {
  28     char buf[32];
  29     string result = "";
  30     for( size_t i = 0; i < nelems; i++ )
  31     {
  32         sprintf(buf, "%d", v[i]);
  33         result += string(buf);
  34         if( i < nelems - 1 )
  35             result += sep;
  36     }
  37     return result;
  38 }
  39
  40
  41 Size randomSize(RNG& rng, double maxSizeLog)
  42 {
  43     double width_log = rng.uniform(0., maxSizeLog);
  44     double height_log = rng.uniform(0., maxSizeLog - width_log);
  45     if( (unsigned)rng % 2 != 0 )
  46         std::swap(width_log, height_log);
  47     Size sz;
  48     sz.width = cvRound(exp(width_log));
  49     sz.height = cvRound(exp(height_log));
  50     return sz;
  51 }
  52
  53 void randomSize(RNG& rng, int minDims, int maxDims, double maxSizeLog, vector<int>& sz)
  54 {
  55     int i, dims = rng.uniform(minDims, maxDims+1);
  56     sz.resize(dims);
  57     for( i = 0; i < dims; i++ )
  58     {
  59         double v = rng.uniform(0., maxSizeLog);
  60         maxSizeLog -= v;
  61         sz[i] = cvRound(exp(v));
  62     }
  63     for( i = 0; i < dims; i++ )
  64     {
  65         int j = rng.uniform(0, dims);
  66         int k = rng.uniform(0, dims);
  67         std::swap(sz[j], sz[k]);
  68     }
  69 }
  70
  71 int randomType(RNG& rng, int typeMask, int minChannels, int maxChannels)
  72 {
  73     int channels = rng.uniform(minChannels, maxChannels+1);
  74     int depth = 0;
  75     CV_Assert((typeMask & _OutputArray::DEPTH_MASK_ALL) != 0);
  76     for(;;)
  77     {
  78         depth = rng.uniform(CV_8U, CV_64F+1);
  79         if( ((1 << depth) & typeMask) != 0 )
  80             break;
  81     }
  82     return CV_MAKETYPE(depth, channels);
  83 }
  84
  85 double getMinVal(int depth)
  86 {
  87     depth = CV_MAT_DEPTH(depth);
  88     double val = depth == CV_8U ? 0 : depth == CV_8S ? SCHAR_MIN : depth == CV_16U ? 0 :
  89     depth == CV_16S ? SHRT_MIN : depth == CV_32S ? INT_MIN :
  90     depth == CV_32F ? -FLT_MAX : depth == CV_64F ? -DBL_MAX : -1;
  91     CV_Assert(val != -1);
  92     return val;
  93 }
  94
  95 double getMaxVal(int depth)
  96 {
  97     depth = CV_MAT_DEPTH(depth);
  98     double val = depth == CV_8U ? UCHAR_MAX : depth == CV_8S ? SCHAR_MAX : depth == CV_16U ? USHRT_MAX :
  99     depth == CV_16S ? SHRT_MAX : depth == CV_32S ? INT_MAX :
 100     depth == CV_32F ? FLT_MAX : depth == CV_64F ? DBL_MAX : -1;
 101     CV_Assert(val != -1);
 102     return val;
 103 }
 104
 105 Mat randomMat(RNG& rng, Size size, int type, double minVal, double maxVal, bool useRoi)
 106 {
 107     Size size0 = size;
 108     if( useRoi )
 109     {
 110         size0.width += std::max(rng.uniform(0, 10) - 5, 0);
 111         size0.height += std::max(rng.uniform(0, 10) - 5, 0);
 112     }
 113
 114     Mat m(size0, type);
 115
 116     rng.fill(m, RNG::UNIFORM, minVal, maxVal);
 117     if( size0 == size )
 118         return m;
 119     return m(Rect((size0.width-size.width)/2, (size0.height-size.height)/2, size.width, size.height));
 120 }
 121
 122 Mat randomMat(RNG& rng, const vector<int>& size, int type, double minVal, double maxVal, bool useRoi)
 123 {
 124     int i, dims = (int)size.size();
 125     vector<int> size0(dims);
 126     vector<Range> r(dims);
 127     bool eqsize = true;
 128     for( i = 0; i < dims; i++ )
 129     {
 130         size0[i] = size[i];
 131         r[i] = Range::all();
 132         if( useRoi )
 133         {
 134             size0[i] += std::max(rng.uniform(0, 5) - 2, 0);
 135             r[i] = Range((size0[i] - size[i])/2, (size0[i] - size[i])/2 + size[i]);
 136         }
 137         eqsize = eqsize && size[i] == size0[i];
 138     }
 139
 140     Mat m(dims, &size0[0], type);
 141
 142     rng.fill(m, RNG::UNIFORM, minVal, maxVal);
 143     if( eqsize )
 144         return m;
 145     return m(&r[0]);
 146 }
 147
 148 void add(const Mat& _a, double alpha, const Mat& _b, double beta,
 149         Scalar gamma, Mat& c, int ctype, bool calcAbs)
 150 {
 151     Mat a = _a, b = _b;
 152     if( a.empty() || alpha == 0 )
 153     {
 154         // both alpha and beta can be 0, but at least one of a and b must be non-empty array,
 155         // otherwise we do not know the size of the output (and may be type of the output, when ctype<0)
 156         CV_Assert( !a.empty() || !b.empty() );
 157         if( !b.empty() )
 158         {
 159             a = b;
 160             alpha = beta;
 161             b = Mat();
 162             beta = 0;
 163         }
 164     }
 165     if( b.empty() || beta == 0 )
 166     {
 167         b = Mat();
 168         beta = 0;
 169     }
 170     else
 171         CV_Assert(a.size == b.size);
 172
 173     if( ctype < 0 )
 174         ctype = a.depth();
 175     ctype = CV_MAKETYPE(CV_MAT_DEPTH(ctype), a.channels());
 176     c.create(a.dims, &a.size[0], ctype);
 177     const Mat *arrays[] = {&a, &b, &c, 0};
 178     Mat planes[3], buf[3];
 179
 180     NAryMatIterator it(arrays, planes);
 181     size_t i, nplanes = it.nplanes;
 182     int cn=a.channels();
 183     int total = (int)planes[0].total(), maxsize = std::min(12*12*std::max(12/cn, 1), total);
 184
 185     CV_Assert(planes[0].rows == 1);
 186     buf[0].create(1, maxsize, CV_64FC(cn));
 187     if(!b.empty())
 188         buf[1].create(1, maxsize, CV_64FC(cn));
 189     buf[2].create(1, maxsize, CV_64FC(cn));
 190     scalarToRawData(gamma, buf[2].ptr(), CV_64FC(cn), (int)(maxsize*cn));
 191
 192     for( i = 0; i < nplanes; i++, ++it)
 193     {
 194         for( int j = 0; j < total; j += maxsize )
 195         {
 196             int j2 = std::min(j + maxsize, total);
 197             Mat apart0 = planes[0].colRange(j, j2);
 198             Mat cpart0 = planes[2].colRange(j, j2);
 199             Mat apart = buf[0].colRange(0, j2 - j);
 200
 201             apart0.convertTo(apart, apart.type(), alpha);
 202             size_t k, n = (j2 - j)*cn;
 203             double* aptr = apart.ptr<double>();
 204             const double* gptr = buf[2].ptr<double>();
 205
 206             if( b.empty() )
 207             {
 208                 for( k = 0; k < n; k++ )
 209                     aptr[k] += gptr[k];
 210             }
 211             else
 212             {
 213                 Mat bpart0 = planes[1].colRange((int)j, (int)j2);
 214                 Mat bpart = buf[1].colRange(0, (int)(j2 - j));
 215                 bpart0.convertTo(bpart, bpart.type(), beta);
 216                 const double* bptr = bpart.ptr<double>();
 217
 218                 for( k = 0; k < n; k++ )
 219                     aptr[k] += bptr[k] + gptr[k];
 220             }
 221             if( calcAbs )
 222                 for( k = 0; k < n; k++ )
 223                     aptr[k] = fabs(aptr[k]);
 224             apart.convertTo(cpart0, cpart0.type(), 1, 0);
 225         }
 226     }
 227 }
 228
 229
 230 template<typename _Tp1, typename _Tp2> inline void
 231 convert_(const _Tp1* src, _Tp2* dst, size_t total, double alpha, double beta)
 232 {
 233     size_t i;
 234     if( alpha == 1 && beta == 0 )
 235         for( i = 0; i < total; i++ )
 236             dst[i] = saturate_cast<_Tp2>(src[i]);
 237     else if( beta == 0 )
 238         for( i = 0; i < total; i++ )
 239             dst[i] = saturate_cast<_Tp2>(src[i]*alpha);
 240     else
 241         for( i = 0; i < total; i++ )
 242             dst[i] = saturate_cast<_Tp2>(src[i]*alpha + beta);
 243 }
 244
 245 template<typename _Tp> inline void
 246 convertTo(const _Tp* src, void* dst, int dtype, size_t total, double alpha, double beta)
 247 {
 248     switch( CV_MAT_DEPTH(dtype) )
 249     {
 250     case CV_8U:
 251         convert_(src, (uchar*)dst, total, alpha, beta);
 252         break;
 253     case CV_8S:
 254         convert_(src, (schar*)dst, total, alpha, beta);
 255         break;
 256     case CV_16U:
 257         convert_(src, (ushort*)dst, total, alpha, beta);
 258         break;
 259     case CV_16S:
 260         convert_(src, (short*)dst, total, alpha, beta);
 261         break;
 262     case CV_32S:
 263         convert_(src, (int*)dst, total, alpha, beta);
 264         break;
 265     case CV_32F:
 266         convert_(src, (float*)dst, total, alpha, beta);
 267         break;
 268     case CV_64F:
 269         convert_(src, (double*)dst, total, alpha, beta);
 270         break;
 271     default:
 272         CV_Assert(0);
 273     }
 274 }
 275
 276 void convert(const Mat& src, cv::OutputArray _dst, int dtype, double alpha, double beta)
 277 {
 278     if (dtype < 0) dtype = _dst.depth();
 279
 280     dtype = CV_MAKETYPE(CV_MAT_DEPTH(dtype), src.channels());
 281     _dst.create(src.dims, &src.size[0], dtype);
 282     Mat dst = _dst.getMat();
 283     if( alpha == 0 )
 284     {
 285         set( dst, Scalar::all(beta) );
 286         return;
 287     }
 288     if( dtype == src.type() && alpha == 1 && beta == 0 )
 289     {
 290         copy( src, dst );
 291         return;
 292     }
 293
 294     const Mat *arrays[]={&src, &dst, 0};
 295     Mat planes[2];
 296
 297     NAryMatIterator it(arrays, planes);
 298     size_t total = planes[0].total()*planes[0].channels();
 299     size_t i, nplanes = it.nplanes;
 300
 301     for( i = 0; i < nplanes; i++, ++it)
 302     {
 303         const uchar* sptr = planes[0].ptr();
 304         uchar* dptr = planes[1].ptr();
 305
 306         switch( src.depth() )
 307         {
 308         case CV_8U:
 309             convertTo((const uchar*)sptr, dptr, dtype, total, alpha, beta);
 310             break;
 311         case CV_8S:
 312             convertTo((const schar*)sptr, dptr, dtype, total, alpha, beta);
 313             break;
 314         case CV_16U:
 315             convertTo((const ushort*)sptr, dptr, dtype, total, alpha, beta);
 316             break;
 317         case CV_16S:
 318             convertTo((const short*)sptr, dptr, dtype, total, alpha, beta);
 319             break;
 320         case CV_32S:
 321             convertTo((const int*)sptr, dptr, dtype, total, alpha, beta);
 322             break;
 323         case CV_32F:
 324             convertTo((const float*)sptr, dptr, dtype, total, alpha, beta);
 325             break;
 326         case CV_64F:
 327             convertTo((const double*)sptr, dptr, dtype, total, alpha, beta);
 328             break;
 329         }
 330     }
 331 }
 332
 333
 334 void copy(const Mat& src, Mat& dst, const Mat& mask, bool invertMask)
 335 {
 336     dst.create(src.dims, &src.size[0], src.type());
 337
 338     if(mask.empty())
 339     {
 340         const Mat* arrays[] = {&src, &dst, 0};
 341         Mat planes[2];
 342         NAryMatIterator it(arrays, planes);
 343         size_t i, nplanes = it.nplanes;
 344         size_t planeSize = planes[0].total()*src.elemSize();
 345
 346         for( i = 0; i < nplanes; i++, ++it )
 347             memcpy(planes[1].ptr(), planes[0].ptr(), planeSize);
 348
 349         return;
 350     }
 351
 352     int mcn = mask.channels();
 353     CV_Assert( src.size == mask.size && mask.depth() == CV_8U
 354                && (mcn == 1 || mcn == src.channels()) );
 355
 356     const Mat *arrays[]={&src, &dst, &mask, 0};
 357     Mat planes[3];
 358
 359     NAryMatIterator it(arrays, planes);
 360     size_t j, k, elemSize = src.elemSize(), maskElemSize = mask.elemSize(), total = planes[0].total();
 361     size_t i, nplanes = it.nplanes;
 362     size_t elemSize1 = src.elemSize1();
 363
 364     for( i = 0; i < nplanes; i++, ++it)
 365     {
 366         const uchar* sptr = planes[0].ptr();
 367         uchar* dptr = planes[1].ptr();
 368         const uchar* mptr = planes[2].ptr();
 369         for( j = 0; j < total; j++, sptr += elemSize, dptr += elemSize, mptr += maskElemSize )
 370         {
 371             if( mcn == 1)
 372             {
 373                 if( (mptr[0] != 0) ^ invertMask )
 374                     for( k = 0; k < elemSize; k++ )
 375                         dptr[k] = sptr[k];
 376             }
 377             else
 378             {
 379                 for( int c = 0; c < mcn; c++ )
 380                     if( (mptr[c] != 0) ^ invertMask )
 381                         for( k = 0; k < elemSize1; k++ )
 382                             dptr[k + c * elemSize1] = sptr[k + c * elemSize1];
 383             }
 384         }
 385     }
 386 }
 387
 388
 389 void set(Mat& dst, const Scalar& gamma, const Mat& mask)
 390 {
 391     double buf[12];
 392     scalarToRawData(gamma, &buf, dst.type(), dst.channels());
 393     const uchar* gptr = (const uchar*)&buf[0];
 394
 395     if(mask.empty())
 396     {
 397         const Mat* arrays[] = {&dst, 0};
 398         Mat plane;
 399         NAryMatIterator it(arrays, &plane);
 400         size_t i, nplanes = it.nplanes;
 401         size_t j, k, elemSize = dst.elemSize(), planeSize = plane.total()*elemSize;
 402
 403         for( k = 1; k < elemSize; k++ )
 404             if( gptr[k] != gptr[0] )
 405                 break;
 406         bool uniform = k >= elemSize;
 407
 408         for( i = 0; i < nplanes; i++, ++it )
 409         {
 410             uchar* dptr = plane.ptr();
 411             if( uniform )
 412                 memset( dptr, gptr[0], planeSize );
 413             else if( i == 0 )
 414             {
 415                 for( j = 0; j < planeSize; j += elemSize, dptr += elemSize )
 416                     for( k = 0; k < elemSize; k++ )
 417                         dptr[k] = gptr[k];
 418             }
 419             else
 420                 memcpy(dptr, dst.ptr(), planeSize);
 421         }
 422         return;
 423     }
 424
 425     int cn = dst.channels(), mcn = mask.channels();
 426     CV_Assert( dst.size == mask.size && (mcn == 1 || mcn == cn) );
 427
 428     const Mat *arrays[]={&dst, &mask, 0};
 429     Mat planes[2];
 430
 431     NAryMatIterator it(arrays, planes);
 432     size_t j, k, elemSize = dst.elemSize(), maskElemSize = mask.elemSize(), total = planes[0].total();
 433     size_t i, nplanes = it.nplanes;
 434     size_t elemSize1 = dst.elemSize1();
 435
 436     for( i = 0; i < nplanes; i++, ++it)
 437     {
 438         uchar* dptr = planes[0].ptr();
 439         const uchar* mptr = planes[1].ptr();
 440
 441         for( j = 0; j < total; j++, dptr += elemSize, mptr += maskElemSize )
 442         {
 443             if( mcn == 1)
 444             {
 445                 if( mptr[0] )
 446                     for( k = 0; k < elemSize; k++ )
 447                         dptr[k] = gptr[k];
 448             }
 449             else
 450             {
 451                 for( int c = 0; c < mcn; c++ )
 452                     if( mptr[c] )
 453                         for( k = 0; k < elemSize1; k++ )
 454                             dptr[k + c * elemSize1] = gptr[k + c * elemSize1];
 455             }
 456         }
 457     }
 458 }
 459
 460
 461 void insert(const Mat& src, Mat& dst, int coi)
 462 {
 463     CV_Assert( dst.size == src.size && src.depth() == dst.depth() &&
 464               0 <= coi && coi < dst.channels() );
 465
 466     const Mat* arrays[] = {&src, &dst, 0};
 467     Mat planes[2];
 468     NAryMatIterator it(arrays, planes);
 469     size_t i, nplanes = it.nplanes;
 470     size_t j, k, size0 = src.elemSize(), size1 = dst.elemSize(), total = planes[0].total();
 471
 472     for( i = 0; i < nplanes; i++, ++it )
 473     {
 474         const uchar* sptr = planes[0].ptr();
 475         uchar* dptr = planes[1].ptr() + coi*size0;
 476
 477         for( j = 0; j < total; j++, sptr += size0, dptr += size1 )
 478         {
 479             for( k = 0; k < size0; k++ )
 480                 dptr[k] = sptr[k];
 481         }
 482     }
 483 }
 484
 485
 486 void extract(const Mat& src, Mat& dst, int coi)
 487 {
 488     dst.create( src.dims, &src.size[0], src.depth() );
 489     CV_Assert( 0 <= coi && coi < src.channels() );
 490
 491     const Mat* arrays[] = {&src, &dst, 0};
 492     Mat planes[2];
 493     NAryMatIterator it(arrays, planes);
 494     size_t i, nplanes = it.nplanes;
 495     size_t j, k, size0 = src.elemSize(), size1 = dst.elemSize(), total = planes[0].total();
 496
 497     for( i = 0; i < nplanes; i++, ++it )
 498     {
 499         const uchar* sptr = planes[0].ptr() + coi*size1;
 500         uchar* dptr = planes[1].ptr();
 501
 502         for( j = 0; j < total; j++, sptr += size0, dptr += size1 )
 503         {
 504             for( k = 0; k < size1; k++ )
 505                 dptr[k] = sptr[k];
 506         }
 507     }
 508 }
 509
 510
 511 void transpose(const Mat& src, Mat& dst)
 512 {
 513     CV_Assert(src.data != dst.data && "Inplace is not support in cvtest::transpose");
 514     CV_Assert(src.dims == 2);
 515     dst.create(src.cols, src.rows, src.type());
 516     int i, j, k, esz = (int)src.elemSize();
 517
 518     for( i = 0; i < dst.rows; i++ )
 519     {
 520         const uchar* sptr = src.ptr(0) + i*esz;
 521         uchar* dptr = dst.ptr(i);
 522
 523         for( j = 0; j < dst.cols; j++, sptr += src.step[0], dptr += esz )
 524         {
 525             for( k = 0; k < esz; k++ )
 526                 dptr[k] = sptr[k];
 527         }
 528     }
 529 }
 530
 531
 532 template<typename _Tp> static void
 533 randUniInt_(RNG& rng, _Tp* data, size_t total, int cn, const Scalar& scale, const Scalar& delta)
 534 {
 535     for( size_t i = 0; i < total; i += cn )
 536         for( int k = 0; k < cn; k++ )
 537         {
 538             int val = cvFloor( randInt(rng)*scale[k] + delta[k] );
 539             data[i + k] = saturate_cast<_Tp>(val);
 540         }
 541 }
 542
 543
 544 template<typename _Tp> static void
 545 randUniFlt_(RNG& rng, _Tp* data, size_t total, int cn, const Scalar& scale, const Scalar& delta)
 546 {
 547     for( size_t i = 0; i < total; i += cn )
 548         for( int k = 0; k < cn; k++ )
 549         {
 550             double val = randReal(rng)*scale[k] + delta[k];
 551             data[i + k] = saturate_cast<_Tp>(val);
 552         }
 553 }
 554
 555
 556 void randUni( RNG& rng, Mat& a, const Scalar& param0, const Scalar& param1 )
 557 {
 558     Scalar scale = param0;
 559     Scalar delta = param1;
 560     double C = a.depth() < CV_32F ? 1./(65536.*65536.) : 1.;
 561
 562     for( int k = 0; k < 4; k++ )
 563     {
 564         double s = scale.val[k] - delta.val[k];
 565         if( s >= 0 )
 566             scale.val[k] = s;
 567         else
 568         {
 569             delta.val[k] = scale.val[k];
 570             scale.val[k] = -s;
 571         }
 572         scale.val[k] *= C;
 573     }
 574
 575     const Mat *arrays[]={&a, 0};
 576     Mat plane;
 577
 578     NAryMatIterator it(arrays, &plane);
 579     size_t i, nplanes = it.nplanes;
 580     int depth = a.depth(), cn = a.channels();
 581     size_t total = plane.total()*cn;
 582
 583     for( i = 0; i < nplanes; i++, ++it )
 584     {
 585         switch( depth )
 586         {
 587         case CV_8U:
 588             randUniInt_(rng, plane.ptr<uchar>(), total, cn, scale, delta);
 589             break;
 590         case CV_8S:
 591             randUniInt_(rng, plane.ptr<schar>(), total, cn, scale, delta);
 592             break;
 593         case CV_16U:
 594             randUniInt_(rng, plane.ptr<ushort>(), total, cn, scale, delta);
 595             break;
 596         case CV_16S:
 597             randUniInt_(rng, plane.ptr<short>(), total, cn, scale, delta);
 598             break;
 599         case CV_32S:
 600             randUniInt_(rng, plane.ptr<int>(), total, cn, scale, delta);
 601             break;
 602         case CV_32F:
 603             randUniFlt_(rng, plane.ptr<float>(), total, cn, scale, delta);
 604             break;
 605         case CV_64F:
 606             randUniFlt_(rng, plane.ptr<double>(), total, cn, scale, delta);
 607             break;
 608         default:
 609             CV_Assert(0);
 610         }
 611     }
 612 }
 613
 614
 615 template<typename _Tp> static void
 616 erode_(const Mat& src, Mat& dst, const vector<int>& ofsvec)
 617 {
 618     int width = dst.cols*src.channels(), n = (int)ofsvec.size();
 619     const int* ofs = &ofsvec[0];
 620
 621     for( int y = 0; y < dst.rows; y++ )
 622     {
 623         const _Tp* sptr = src.ptr<_Tp>(y);
 624         _Tp* dptr = dst.ptr<_Tp>(y);
 625
 626         for( int x = 0; x < width; x++ )
 627         {
 628             _Tp result = sptr[x + ofs[0]];
 629             for( int i = 1; i < n; i++ )
 630                 result = std::min(result, sptr[x + ofs[i]]);
 631             dptr[x] = result;
 632         }
 633     }
 634 }
 635
 636
 637 template<typename _Tp> static void
 638 dilate_(const Mat& src, Mat& dst, const vector<int>& ofsvec)
 639 {
 640     int width = dst.cols*src.channels(), n = (int)ofsvec.size();
 641     const int* ofs = &ofsvec[0];
 642
 643     for( int y = 0; y < dst.rows; y++ )
 644     {
 645         const _Tp* sptr = src.ptr<_Tp>(y);
 646         _Tp* dptr = dst.ptr<_Tp>(y);
 647
 648         for( int x = 0; x < width; x++ )
 649         {
 650             _Tp result = sptr[x + ofs[0]];
 651             for( int i = 1; i < n; i++ )
 652                 result = std::max(result, sptr[x + ofs[i]]);
 653             dptr[x] = result;
 654         }
 655     }
 656 }
 657
 658
 659 void erode(const Mat& _src, Mat& dst, const Mat& _kernel, Point anchor,
 660            int borderType, const Scalar& _borderValue)
 661 {
 662     //if( _src.type() == CV_16UC3 && _src.size() == Size(1, 2) )
 663     //    putchar('*');
 664     Mat kernel = _kernel, src;
 665     Scalar borderValue = _borderValue;
 666     if( kernel.empty() )
 667         kernel = Mat::ones(3, 3, CV_8U);
 668     else
 669     {
 670         CV_Assert( kernel.type() == CV_8U );
 671     }
 672     if( anchor == Point(-1,-1) )
 673         anchor = Point(kernel.cols/2, kernel.rows/2);
 674     if( borderType == BORDER_CONSTANT )
 675         borderValue = getMaxVal(src.depth());
 676     copyMakeBorder(_src, src, anchor.y, kernel.rows - anchor.y - 1,
 677                    anchor.x, kernel.cols - anchor.x - 1,
 678                    borderType, borderValue);
 679     dst.create( _src.size(), src.type() );
 680
 681     vector<int> ofs;
 682     int step = (int)(src.step/src.elemSize1()), cn = src.channels();
 683     for( int i = 0; i < kernel.rows; i++ )
 684         for( int j = 0; j < kernel.cols; j++ )
 685             if( kernel.at<uchar>(i, j) != 0 )
 686                 ofs.push_back(i*step + j*cn);
 687     if( ofs.empty() )
 688         ofs.push_back(anchor.y*step + anchor.x*cn);
 689
 690     switch( src.depth() )
 691     {
 692     case CV_8U:
 693         erode_<uchar>(src, dst, ofs);
 694         break;
 695     case CV_8S:
 696         erode_<schar>(src, dst, ofs);
 697         break;
 698     case CV_16U:
 699         erode_<ushort>(src, dst, ofs);
 700         break;
 701     case CV_16S:
 702         erode_<short>(src, dst, ofs);
 703         break;
 704     case CV_32S:
 705         erode_<int>(src, dst, ofs);
 706         break;
 707     case CV_32F:
 708         erode_<float>(src, dst, ofs);
 709         break;
 710     case CV_64F:
 711         erode_<double>(src, dst, ofs);
 712         break;
 713     default:
 714         CV_Assert(0);
 715     }
 716 }
 717
 718 void dilate(const Mat& _src, Mat& dst, const Mat& _kernel, Point anchor,
 719             int borderType, const Scalar& _borderValue)
 720 {
 721     Mat kernel = _kernel, src;
 722     Scalar borderValue = _borderValue;
 723     if( kernel.empty() )
 724         kernel = Mat::ones(3, 3, CV_8U);
 725     else
 726     {
 727         CV_Assert( kernel.type() == CV_8U );
 728     }
 729     if( anchor == Point(-1,-1) )
 730         anchor = Point(kernel.cols/2, kernel.rows/2);
 731     if( borderType == BORDER_CONSTANT )
 732         borderValue = getMinVal(src.depth());
 733     copyMakeBorder(_src, src, anchor.y, kernel.rows - anchor.y - 1,
 734                    anchor.x, kernel.cols - anchor.x - 1,
 735                    borderType, borderValue);
 736     dst.create( _src.size(), src.type() );
 737
 738     vector<int> ofs;
 739     int step = (int)(src.step/src.elemSize1()), cn = src.channels();
 740     for( int i = 0; i < kernel.rows; i++ )
 741         for( int j = 0; j < kernel.cols; j++ )
 742             if( kernel.at<uchar>(i, j) != 0 )
 743                 ofs.push_back(i*step + j*cn);
 744     if( ofs.empty() )
 745         ofs.push_back(anchor.y*step + anchor.x*cn);
 746
 747     switch( src.depth() )
 748     {
 749     case CV_8U:
 750         dilate_<uchar>(src, dst, ofs);
 751         break;
 752     case CV_8S:
 753         dilate_<schar>(src, dst, ofs);
 754         break;
 755     case CV_16U:
 756         dilate_<ushort>(src, dst, ofs);
 757         break;
 758     case CV_16S:
 759         dilate_<short>(src, dst, ofs);
 760         break;
 761     case CV_32S:
 762         dilate_<int>(src, dst, ofs);
 763         break;
 764     case CV_32F:
 765         dilate_<float>(src, dst, ofs);
 766         break;
 767     case CV_64F:
 768         dilate_<double>(src, dst, ofs);
 769         break;
 770     default:
 771         CV_Assert(0);
 772     }
 773 }
 774
 775
 776 template<typename _Tp> static void
 777 filter2D_(const Mat& src, Mat& dst, const vector<int>& ofsvec, const vector<double>& coeffvec)
 778 {
 779     const int* ofs = &ofsvec[0];
 780     const double* coeff = &coeffvec[0];
 781     int width = dst.cols*dst.channels(), ncoeffs = (int)ofsvec.size();
 782
 783     for( int y = 0; y < dst.rows; y++ )
 784     {
 785         const _Tp* sptr = src.ptr<_Tp>(y);
 786         double* dptr = dst.ptr<double>(y);
 787
 788         for( int x = 0; x < width; x++ )
 789         {
 790             double s = 0;
 791             for( int i = 0; i < ncoeffs; i++ )
 792                 s += sptr[x + ofs[i]]*coeff[i];
 793             dptr[x] = s;
 794         }
 795     }
 796 }
 797
 798
 799 void filter2D(const Mat& _src, Mat& dst, int ddepth, const Mat& kernel,
 800               Point anchor, double delta, int borderType, const Scalar& _borderValue)
 801 {
 802     Mat src, _dst;
 803     Scalar borderValue = _borderValue;
 804     CV_Assert( kernel.type() == CV_32F || kernel.type() == CV_64F );
 805     if( anchor == Point(-1,-1) )
 806         anchor = Point(kernel.cols/2, kernel.rows/2);
 807     if( borderType == BORDER_CONSTANT )
 808         borderValue = getMinVal(src.depth());
 809     copyMakeBorder(_src, src, anchor.y, kernel.rows - anchor.y - 1,
 810                    anchor.x, kernel.cols - anchor.x - 1,
 811                    borderType, borderValue);
 812     _dst.create( _src.size(), CV_MAKETYPE(CV_64F, src.channels()) );
 813
 814     vector<int> ofs;
 815     vector<double> coeff(kernel.rows*kernel.cols);
 816     Mat cmat(kernel.rows, kernel.cols, CV_64F, &coeff[0]);
 817     convert(kernel, cmat, cmat.type());
 818
 819     int step = (int)(src.step/src.elemSize1()), cn = src.channels();
 820     for( int i = 0; i < kernel.rows; i++ )
 821         for( int j = 0; j < kernel.cols; j++ )
 822                 ofs.push_back(i*step + j*cn);
 823
 824     switch( src.depth() )
 825     {
 826     case CV_8U:
 827         filter2D_<uchar>(src, _dst, ofs, coeff);
 828         break;
 829     case CV_8S:
 830         filter2D_<schar>(src, _dst, ofs, coeff);
 831         break;
 832     case CV_16U:
 833         filter2D_<ushort>(src, _dst, ofs, coeff);
 834         break;
 835     case CV_16S:
 836         filter2D_<short>(src, _dst, ofs, coeff);
 837         break;
 838     case CV_32S:
 839         filter2D_<int>(src, _dst, ofs, coeff);
 840         break;
 841     case CV_32F:
 842         filter2D_<float>(src, _dst, ofs, coeff);
 843         break;
 844     case CV_64F:
 845         filter2D_<double>(src, _dst, ofs, coeff);
 846         break;
 847     default:
 848         CV_Assert(0);
 849     }
 850
 851     convert(_dst, dst, ddepth, 1, delta);
 852 }
 853
 854
 855 static int borderInterpolate( int p, int len, int borderType )
 856 {
 857     if( (unsigned)p < (unsigned)len )
 858         ;
 859     else if( borderType == BORDER_REPLICATE )
 860         p = p < 0 ? 0 : len - 1;
 861     else if( borderType == BORDER_REFLECT || borderType == BORDER_REFLECT_101 )
 862     {
 863         int delta = borderType == BORDER_REFLECT_101;
 864         if( len == 1 )
 865             return 0;
 866         do
 867         {
 868             if( p < 0 )
 869                 p = -p - 1 + delta;
 870             else
 871                 p = len - 1 - (p - len) - delta;
 872         }
 873         while( (unsigned)p >= (unsigned)len );
 874     }
 875     else if( borderType == BORDER_WRAP )
 876     {
 877         if( p < 0 )
 878             p -= ((p-len+1)/len)*len;
 879         if( p >= len )
 880             p %= len;
 881     }
 882     else if( borderType == BORDER_CONSTANT )
 883         p = -1;
 884     else
 885         CV_Error( Error::StsBadArg, "Unknown/unsupported border type" );
 886     return p;
 887 }
 888
 889
 890 void copyMakeBorder(const Mat& src, Mat& dst, int top, int bottom, int left, int right,
 891                     int borderType, const Scalar& borderValue)
 892 {
 893     dst.create(src.rows + top + bottom, src.cols + left + right, src.type());
 894     int i, j, k, esz = (int)src.elemSize();
 895     int width = src.cols*esz, width1 = dst.cols*esz;
 896
 897     if( borderType == BORDER_CONSTANT )
 898     {
 899         vector<uchar> valvec((src.cols + left + right)*esz);
 900         uchar* val = &valvec[0];
 901         scalarToRawData(borderValue, val, src.type(), (src.cols + left + right)*src.channels());
 902
 903         left *= esz;
 904         right *= esz;
 905         for( i = 0; i < src.rows; i++ )
 906         {
 907             const uchar* sptr = src.ptr(i);
 908             uchar* dptr = dst.ptr(i + top) + left;
 909             for( j = 0; j < left; j++ )
 910                 dptr[j - left] = val[j];
 911             if( dptr != sptr )
 912                 for( j = 0; j < width; j++ )
 913                     dptr[j] = sptr[j];
 914             for( j = 0; j < right; j++ )
 915                 dptr[j + width] = val[j];
 916         }
 917
 918         for( i = 0; i < top; i++ )
 919         {
 920             uchar* dptr = dst.ptr(i);
 921             for( j = 0; j < width1; j++ )
 922                 dptr[j] = val[j];
 923         }
 924
 925         for( i = 0; i < bottom; i++ )
 926         {
 927             uchar* dptr = dst.ptr(i + top + src.rows);
 928             for( j = 0; j < width1; j++ )
 929                 dptr[j] = val[j];
 930         }
 931     }
 932     else
 933     {
 934         vector<int> tabvec((left + right)*esz + 1);
 935         int* ltab = &tabvec[0];
 936         int* rtab = &tabvec[left*esz];
 937         for( i = 0; i < left; i++ )
 938         {
 939             j = borderInterpolate(i - left, src.cols, borderType)*esz;
 940             for( k = 0; k < esz; k++ )
 941                 ltab[i*esz + k] = j + k;
 942         }
 943         for( i = 0; i < right; i++ )
 944         {
 945             j = borderInterpolate(src.cols + i, src.cols, borderType)*esz;
 946             for( k = 0; k < esz; k++ )
 947                 rtab[i*esz + k] = j + k;
 948         }
 949
 950         left *= esz;
 951         right *= esz;
 952         for( i = 0; i < src.rows; i++ )
 953         {
 954             const uchar* sptr = src.ptr(i);
 955             uchar* dptr = dst.ptr(i + top);
 956
 957             for( j = 0; j < left; j++ )
 958                 dptr[j] = sptr[ltab[j]];
 959             if( dptr + left != sptr )
 960             {
 961                 for( j = 0; j < width; j++ )
 962                     dptr[j + left] = sptr[j];
 963             }
 964             for( j = 0; j < right; j++ )
 965                 dptr[j + left + width] = sptr[rtab[j]];
 966         }
 967
 968         for( i = 0; i < top; i++ )
 969         {
 970             j = borderInterpolate(i - top, src.rows, borderType);
 971             const uchar* sptr = dst.ptr(j + top);
 972             uchar* dptr = dst.ptr(i);
 973
 974             for( k = 0; k < width1; k++ )
 975                 dptr[k] = sptr[k];
 976         }
 977
 978         for( i = 0; i < bottom; i++ )
 979         {
 980             j = borderInterpolate(i + src.rows, src.rows, borderType);
 981             const uchar* sptr = dst.ptr(j + top);
 982             uchar* dptr = dst.ptr(i + top + src.rows);
 983
 984             for( k = 0; k < width1; k++ )
 985                 dptr[k] = sptr[k];
 986         }
 987     }
 988 }
 989
 990
 991 template<typename _Tp> static void
 992 minMaxLoc_(const _Tp* src, size_t total, size_t startidx,
 993            double* _minval, double* _maxval,
 994            size_t* _minpos, size_t* _maxpos,
 995            const uchar* mask)
 996 {
 997     _Tp maxval = saturate_cast<_Tp>(*_maxval), minval = saturate_cast<_Tp>(*_minval);
 998     size_t minpos = *_minpos, maxpos = *_maxpos;
 999
1000     if( !mask )
1001     {
1002         for( size_t i = 0; i < total; i++ )
1003         {
1004             _Tp val = src[i];
1005             if( minval > val || !minpos )
1006             {
1007                 minval = val;
1008                 minpos = startidx + i;
1009             }
1010             if( maxval < val || !maxpos )
1011             {
1012                 maxval = val;
1013                 maxpos = startidx + i;
1014             }
1015         }
1016     }
1017     else
1018     {
1019         for( size_t i = 0; i < total; i++ )
1020         {
1021             _Tp val = src[i];
1022             if( (minval > val || !minpos) && mask[i] )
1023             {
1024                 minval = val;
1025                 minpos = startidx + i;
1026             }
1027             if( (maxval < val || !maxpos) && mask[i] )
1028             {
1029                 maxval = val;
1030                 maxpos = startidx + i;
1031             }
1032         }
1033     }
1034
1035     *_maxval = maxval;
1036     *_minval = minval;
1037     *_maxpos = maxpos;
1038     *_minpos = minpos;
1039 }
1040
1041
1042 static void setpos( const Mat& mtx, vector<int>& pos, size_t idx )
1043 {
1044     pos.resize(mtx.dims);
1045     if( idx > 0 )
1046     {
1047         idx--;
1048         for( int i = mtx.dims-1; i >= 0; i-- )
1049         {
1050             int sz = mtx.size[i]*(i == mtx.dims-1 ? mtx.channels() : 1);
1051             pos[i] = (int)(idx % sz);
1052             idx /= sz;
1053         }
1054     }
1055     else
1056     {
1057         for( int i = mtx.dims-1; i >= 0; i-- )
1058             pos[i] = -1;
1059     }
1060 }
1061
1062 void minMaxLoc(const Mat& src, double* _minval, double* _maxval,
1063                vector<int>* _minloc, vector<int>* _maxloc,
1064                const Mat& mask)
1065 {
1066     CV_Assert( src.channels() == 1 );
1067     const Mat *arrays[]={&src, &mask, 0};
1068     Mat planes[2];
1069
1070     NAryMatIterator it(arrays, planes);
1071     size_t startidx = 1, total = planes[0].total();
1072     size_t i, nplanes = it.nplanes;
1073     int depth = src.depth();
1074     double minval = 0;
1075     double maxval = 0;
1076     size_t maxidx = 0, minidx = 0;
1077
1078     for( i = 0; i < nplanes; i++, ++it, startidx += total )
1079     {
1080         const uchar* sptr = planes[0].ptr();
1081         const uchar* mptr = planes[1].ptr();
1082
1083         switch( depth )
1084         {
1085         case CV_8U:
1086             minMaxLoc_((const uchar*)sptr, total, startidx,
1087                        &minval, &maxval, &minidx, &maxidx, mptr);
1088             break;
1089         case CV_8S:
1090             minMaxLoc_((const schar*)sptr, total, startidx,
1091                        &minval, &maxval, &minidx, &maxidx, mptr);
1092             break;
1093         case CV_16U:
1094             minMaxLoc_((const ushort*)sptr, total, startidx,
1095                        &minval, &maxval, &minidx, &maxidx, mptr);
1096             break;
1097         case CV_16S:
1098             minMaxLoc_((const short*)sptr, total, startidx,
1099                        &minval, &maxval, &minidx, &maxidx, mptr);
1100             break;
1101         case CV_32S:
1102             minMaxLoc_((const int*)sptr, total, startidx,
1103                        &minval, &maxval, &minidx, &maxidx, mptr);
1104             break;
1105         case CV_32F:
1106             minMaxLoc_((const float*)sptr, total, startidx,
1107                        &minval, &maxval, &minidx, &maxidx, mptr);
1108             break;
1109         case CV_64F:
1110             minMaxLoc_((const double*)sptr, total, startidx,
1111                        &minval, &maxval, &minidx, &maxidx, mptr);
1112             break;
1113         default:
1114             CV_Assert(0);
1115         }
1116     }
1117
1118     if( _maxval )
1119         *_maxval = maxval;
1120     if( _minval )
1121         *_minval = minval;
1122     if( _maxloc )
1123         setpos( src, *_maxloc, maxidx );
1124     if( _minloc )
1125         setpos( src, *_minloc, minidx );
1126 }
1127
1128
1129 static int
1130 normHamming(const uchar* src, size_t total, int cellSize)
1131 {
1132     int result = 0;
1133     int mask = cellSize == 1 ? 1 : cellSize == 2 ? 3 : cellSize == 4 ? 15 : -1;
1134     CV_Assert( mask >= 0 );
1135
1136     for( size_t i = 0; i < total; i++ )
1137     {
1138         unsigned a = src[i];
1139         for( ; a != 0; a >>= cellSize )
1140             result += (a & mask) != 0;
1141     }
1142     return result;
1143 }
1144
1145
1146 template<typename _Tp> static double
1147 norm_(const _Tp* src, size_t total, int cn, int normType, double startval, const uchar* mask)
1148 {
1149     size_t i;
1150     double result = startval;
1151     if( !mask )
1152         total *= cn;
1153
1154     if( normType == NORM_INF )
1155     {
1156         if( !mask )
1157             for( i = 0; i < total; i++ )
1158                 result = std::max(result, (double)std::abs(0+src[i]));// trick with 0 used to quiet gcc warning
1159         else
1160             for( int c = 0; c < cn; c++ )
1161             {
1162                 for( i = 0; i < total; i++ )
1163                     if( mask[i] )
1164                         result = std::max(result, (double)std::abs(0+src[i*cn + c]));
1165             }
1166     }
1167     else if( normType == NORM_L1 )
1168     {
1169         if( !mask )
1170             for( i = 0; i < total; i++ )
1171                 result += std::abs(0+src[i]);
1172         else
1173             for( int c = 0; c < cn; c++ )
1174             {
1175                 for( i = 0; i < total; i++ )
1176                     if( mask[i] )
1177                         result += std::abs(0+src[i*cn + c]);
1178             }
1179     }
1180     else
1181     {
1182         if( !mask )
1183             for( i = 0; i < total; i++ )
1184             {
1185                 double v = src[i];
1186                 result += v*v;
1187             }
1188         else
1189             for( int c = 0; c < cn; c++ )
1190             {
1191                 for( i = 0; i < total; i++ )
1192                     if( mask[i] )
1193                     {
1194                         double v = src[i*cn + c];
1195                         result += v*v;
1196                     }
1197             }
1198     }
1199     return result;
1200 }
1201
1202
1203 template<typename _Tp> static double
1204 norm_(const _Tp* src1, const _Tp* src2, size_t total, int cn, int normType, double startval, const uchar* mask)
1205 {
1206     size_t i;
1207     double result = startval;
1208     if( !mask )
1209         total *= cn;
1210
1211     if( normType == NORM_INF )
1212     {
1213         if( !mask )
1214             for( i = 0; i < total; i++ )
1215                 result = std::max(result, (double)std::abs(src1[i] - src2[i]));
1216         else
1217             for( int c = 0; c < cn; c++ )
1218             {
1219                 for( i = 0; i < total; i++ )
1220                     if( mask[i] )
1221                         result = std::max(result, (double)std::abs(src1[i*cn + c] - src2[i*cn + c]));
1222             }
1223     }
1224     else if( normType == NORM_L1 )
1225     {
1226         if( !mask )
1227             for( i = 0; i < total; i++ )
1228                 result += std::abs(src1[i] - src2[i]);
1229         else
1230             for( int c = 0; c < cn; c++ )
1231             {
1232                 for( i = 0; i < total; i++ )
1233                     if( mask[i] )
1234                         result += std::abs(src1[i*cn + c] - src2[i*cn + c]);
1235             }
1236     }
1237     else
1238     {
1239         if( !mask )
1240             for( i = 0; i < total; i++ )
1241             {
1242                 double v = src1[i] - src2[i];
1243                 result += v*v;
1244             }
1245         else
1246             for( int c = 0; c < cn; c++ )
1247             {
1248                 for( i = 0; i < total; i++ )
1249                     if( mask[i] )
1250                     {
1251                         double v = src1[i*cn + c] - src2[i*cn + c];
1252                         result += v*v;
1253                     }
1254             }
1255     }
1256     return result;
1257 }
1258
1259
1260 double norm(InputArray _src, int normType, InputArray _mask)
1261 {
1262     Mat src = _src.getMat(), mask = _mask.getMat();
1263     if( normType == NORM_HAMMING || normType == NORM_HAMMING2 )
1264     {
1265         if( !mask.empty() )
1266         {
1267             Mat temp;
1268             bitwise_and(src, mask, temp);
1269             return cvtest::norm(temp, normType, Mat());
1270         }
1271
1272         CV_Assert( src.depth() == CV_8U );
1273
1274         const Mat *arrays[]={&src, 0};
1275         Mat planes[1];
1276
1277         NAryMatIterator it(arrays, planes);
1278         size_t total = planes[0].total();
1279         size_t i, nplanes = it.nplanes;
1280         double result = 0;
1281         int cellSize = normType == NORM_HAMMING ? 1 : 2;
1282
1283         for( i = 0; i < nplanes; i++, ++it )
1284             result += normHamming(planes[0].ptr(), total, cellSize);
1285         return result;
1286     }
1287     int normType0 = normType;
1288     normType = normType == NORM_L2SQR ? NORM_L2 : normType;
1289
1290     CV_Assert( mask.empty() || (src.size == mask.size && mask.type() == CV_8U) );
1291     CV_Assert( normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 );
1292
1293     const Mat *arrays[]={&src, &mask, 0};
1294     Mat planes[2];
1295
1296     NAryMatIterator it(arrays, planes);
1297     size_t total = planes[0].total();
1298     size_t i, nplanes = it.nplanes;
1299     int depth = src.depth(), cn = planes[0].channels();
1300     double result = 0;
1301
1302     for( i = 0; i < nplanes; i++, ++it )
1303     {
1304         const uchar* sptr = planes[0].ptr();
1305         const uchar* mptr = planes[1].ptr();
1306
1307         switch( depth )
1308         {
1309         case CV_8U:
1310             result = norm_((const uchar*)sptr, total, cn, normType, result, mptr);
1311             break;
1312         case CV_8S:
1313             result = norm_((const schar*)sptr, total, cn, normType, result, mptr);
1314             break;
1315         case CV_16U:
1316             result = norm_((const ushort*)sptr, total, cn, normType, result, mptr);
1317             break;
1318         case CV_16S:
1319             result = norm_((const short*)sptr, total, cn, normType, result, mptr);
1320             break;
1321         case CV_32S:
1322             result = norm_((const int*)sptr, total, cn, normType, result, mptr);
1323             break;
1324         case CV_32F:
1325             result = norm_((const float*)sptr, total, cn, normType, result, mptr);
1326             break;
1327         case CV_64F:
1328             result = norm_((const double*)sptr, total, cn, normType, result, mptr);
1329             break;
1330         default:
1331             CV_Error(Error::StsUnsupportedFormat, "");
1332         };
1333     }
1334     if( normType0 == NORM_L2 )
1335         result = sqrt(result);
1336     return result;
1337 }
1338
1339
1340 double norm(InputArray _src1, InputArray _src2, int normType, InputArray _mask)
1341 {
1342     Mat src1 = _src1.getMat(), src2 = _src2.getMat(), mask = _mask.getMat();
1343     bool isRelative = (normType & NORM_RELATIVE) != 0;
1344     normType &= ~NORM_RELATIVE;
1345
1346     if( normType == NORM_HAMMING || normType == NORM_HAMMING2 )
1347     {
1348         Mat temp;
1349         bitwise_xor(src1, src2, temp);
1350         if( !mask.empty() )
1351             bitwise_and(temp, mask, temp);
1352
1353         CV_Assert( temp.depth() == CV_8U );
1354
1355         const Mat *arrays[]={&temp, 0};
1356         Mat planes[1];
1357
1358         NAryMatIterator it(arrays, planes);
1359         size_t total = planes[0].total();
1360         size_t i, nplanes = it.nplanes;
1361         double result = 0;
1362         int cellSize = normType == NORM_HAMMING ? 1 : 2;
1363
1364         for( i = 0; i < nplanes; i++, ++it )
1365             result += normHamming(planes[0].ptr(), total, cellSize);
1366         return result;
1367     }
1368     int normType0 = normType;
1369     normType = normType == NORM_L2SQR ? NORM_L2 : normType;
1370
1371     CV_Assert( src1.type() == src2.type() && src1.size == src2.size );
1372     CV_Assert( mask.empty() || (src1.size == mask.size && mask.type() == CV_8U) );
1373     CV_Assert( normType == NORM_INF || normType == NORM_L1 || normType == NORM_L2 );
1374     const Mat *arrays[]={&src1, &src2, &mask, 0};
1375     Mat planes[3];
1376
1377     NAryMatIterator it(arrays, planes);
1378     size_t total = planes[0].total();
1379     size_t i, nplanes = it.nplanes;
1380     int depth = src1.depth(), cn = planes[0].channels();
1381     double result = 0;
1382
1383     for( i = 0; i < nplanes; i++, ++it )
1384     {
1385         const uchar* sptr1 = planes[0].ptr();
1386         const uchar* sptr2 = planes[1].ptr();
1387         const uchar* mptr = planes[2].ptr();
1388
1389         switch( depth )
1390         {
1391         case CV_8U:
1392             result = norm_((const uchar*)sptr1, (const uchar*)sptr2, total, cn, normType, result, mptr);
1393             break;
1394         case CV_8S:
1395             result = norm_((const schar*)sptr1, (const schar*)sptr2, total, cn, normType, result, mptr);
1396             break;
1397         case CV_16U:
1398             result = norm_((const ushort*)sptr1, (const ushort*)sptr2, total, cn, normType, result, mptr);
1399             break;
1400         case CV_16S:
1401             result = norm_((const short*)sptr1, (const short*)sptr2, total, cn, normType, result, mptr);
1402             break;
1403         case CV_32S:
1404             result = norm_((const int*)sptr1, (const int*)sptr2, total, cn, normType, result, mptr);
1405             break;
1406         case CV_32F:
1407             result = norm_((const float*)sptr1, (const float*)sptr2, total, cn, normType, result, mptr);
1408             break;
1409         case CV_64F:
1410             result = norm_((const double*)sptr1, (const double*)sptr2, total, cn, normType, result, mptr);
1411             break;
1412         default:
1413             CV_Error(Error::StsUnsupportedFormat, "");
1414         };
1415     }
1416     if( normType0 == NORM_L2 )
1417         result = sqrt(result);
1418     return isRelative ? result / (cvtest::norm(src2, normType) + DBL_EPSILON) : result;
1419 }
1420
1421 double PSNR(InputArray _src1, InputArray _src2)
1422 {
1423     CV_Assert( _src1.depth() == CV_8U );
1424     double diff = std::sqrt(cvtest::norm(_src1, _src2, NORM_L2SQR)/(_src1.total()*_src1.channels()));
1425     return 20*log10(255./(diff+DBL_EPSILON));
1426 }
1427
1428 template<typename _Tp> static double
1429 crossCorr_(const _Tp* src1, const _Tp* src2, size_t total)
1430 {
1431     double result = 0;
1432     for( size_t i = 0; i < total; i++ )
1433         result += (double)src1[i]*src2[i];
1434     return result;
1435 }
1436
1437 double crossCorr(const Mat& src1, const Mat& src2)
1438 {
1439     CV_Assert( src1.size == src2.size && src1.type() == src2.type() );
1440     const Mat *arrays[]={&src1, &src2, 0};
1441     Mat planes[2];
1442
1443     NAryMatIterator it(arrays, planes);
1444     size_t total = planes[0].total()*planes[0].channels();
1445     size_t i, nplanes = it.nplanes;
1446     int depth = src1.depth();
1447     double result = 0;
1448
1449     for( i = 0; i < nplanes; i++, ++it )
1450     {
1451         const uchar* sptr1 = planes[0].ptr();
1452         const uchar* sptr2 = planes[1].ptr();
1453
1454         switch( depth )
1455         {
1456         case CV_8U:
1457             result += crossCorr_((const uchar*)sptr1, (const uchar*)sptr2, total);
1458             break;
1459         case CV_8S:
1460             result += crossCorr_((const schar*)sptr1, (const schar*)sptr2, total);
1461             break;
1462         case CV_16U:
1463             result += crossCorr_((const ushort*)sptr1, (const ushort*)sptr2, total);
1464             break;
1465         case CV_16S:
1466             result += crossCorr_((const short*)sptr1, (const short*)sptr2, total);
1467             break;
1468         case CV_32S:
1469             result += crossCorr_((const int*)sptr1, (const int*)sptr2, total);
1470             break;
1471         case CV_32F:
1472             result += crossCorr_((const float*)sptr1, (const float*)sptr2, total);
1473             break;
1474         case CV_64F:
1475             result += crossCorr_((const double*)sptr1, (const double*)sptr2, total);
1476             break;
1477         default:
1478             CV_Error(Error::StsUnsupportedFormat, "");
1479         };
1480     }
1481     return result;
1482 }
1483
1484
1485 static void
1486 logicOp_(const uchar* src1, const uchar* src2, uchar* dst, size_t total, char c)
1487 {
1488     size_t i;
1489     if( c == '&' )
1490         for( i = 0; i < total; i++ )
1491             dst[i] = src1[i] & src2[i];
1492     else if( c == '|' )
1493         for( i = 0; i < total; i++ )
1494             dst[i] = src1[i] | src2[i];
1495     else
1496         for( i = 0; i < total; i++ )
1497             dst[i] = src1[i] ^ src2[i];
1498 }
1499
1500 static void
1501 logicOpS_(const uchar* src, const uchar* scalar, uchar* dst, size_t total, char c)
1502 {
1503     const size_t blockSize = 96;
1504     size_t i, j;
1505     if( c == '&' )
1506         for( i = 0; i < total; i += blockSize, dst += blockSize, src += blockSize )
1507         {
1508             size_t sz = MIN(total - i, blockSize);
1509             for( j = 0; j < sz; j++ )
1510                 dst[j] = src[j] & scalar[j];
1511         }
1512     else if( c == '|' )
1513         for( i = 0; i < total; i += blockSize, dst += blockSize, src += blockSize )
1514         {
1515             size_t sz = MIN(total - i, blockSize);
1516             for( j = 0; j < sz; j++ )
1517                 dst[j] = src[j] | scalar[j];
1518         }
1519     else if( c == '^' )
1520     {
1521         for( i = 0; i < total; i += blockSize, dst += blockSize, src += blockSize )
1522         {
1523             size_t sz = MIN(total - i, blockSize);
1524             for( j = 0; j < sz; j++ )
1525                 dst[j] = src[j] ^ scalar[j];
1526         }
1527     }
1528     else
1529         for( i = 0; i < total; i++ )
1530             dst[i] = ~src[i];
1531 }
1532
1533
1534 void logicOp( const Mat& src1, const Mat& src2, Mat& dst, char op )
1535 {
1536     CV_Assert( op == '&' || op == '|' || op == '^' );
1537     CV_Assert( src1.type() == src2.type() && src1.size == src2.size );
1538     dst.create( src1.dims, &src1.size[0], src1.type() );
1539     const Mat *arrays[]={&src1, &src2, &dst, 0};
1540     Mat planes[3];
1541
1542     NAryMatIterator it(arrays, planes);
1543     size_t total = planes[0].total()*planes[0].elemSize();
1544     size_t i, nplanes = it.nplanes;
1545
1546     for( i = 0; i < nplanes; i++, ++it )
1547     {
1548         const uchar* sptr1 = planes[0].ptr();
1549         const uchar* sptr2 = planes[1].ptr();
1550         uchar* dptr = planes[2].ptr();
1551
1552         logicOp_(sptr1, sptr2, dptr, total, op);
1553     }
1554 }
1555
1556
1557 void logicOp(const Mat& src, const Scalar& s, Mat& dst, char op)
1558 {
1559     CV_Assert( op == '&' || op == '|' || op == '^' || op == '~' );
1560     dst.create( src.dims, &src.size[0], src.type() );
1561     const Mat *arrays[]={&src, &dst, 0};
1562     Mat planes[2];
1563
1564     NAryMatIterator it(arrays, planes);
1565     size_t total = planes[0].total()*planes[0].elemSize();
1566     size_t i, nplanes = it.nplanes;
1567     double buf[12];
1568     scalarToRawData(s, buf, src.type(), (int)(96/planes[0].elemSize1()));
1569
1570     for( i = 0; i < nplanes; i++, ++it )
1571     {
1572         const uchar* sptr = planes[0].ptr();
1573         uchar* dptr = planes[1].ptr();
1574
1575         logicOpS_(sptr, (uchar*)&buf[0], dptr, total, op);
1576     }
1577 }
1578
1579
1580 template<typename _Tp> static void
1581 compare_(const _Tp* src1, const _Tp* src2, uchar* dst, size_t total, int cmpop)
1582 {
1583     size_t i;
1584     switch( cmpop )
1585     {
1586     case CMP_LT:
1587         for( i = 0; i < total; i++ )
1588             dst[i] = src1[i] < src2[i] ? 255 : 0;
1589         break;
1590     case CMP_LE:
1591         for( i = 0; i < total; i++ )
1592             dst[i] = src1[i] <= src2[i] ? 255 : 0;
1593         break;
1594     case CMP_EQ:
1595         for( i = 0; i < total; i++ )
1596             dst[i] = src1[i] == src2[i] ? 255 : 0;
1597         break;
1598     case CMP_NE:
1599         for( i = 0; i < total; i++ )
1600             dst[i] = src1[i] != src2[i] ? 255 : 0;
1601         break;
1602     case CMP_GE:
1603         for( i = 0; i < total; i++ )
1604             dst[i] = src1[i] >= src2[i] ? 255 : 0;
1605         break;
1606     case CMP_GT:
1607         for( i = 0; i < total; i++ )
1608             dst[i] = src1[i] > src2[i] ? 255 : 0;
1609         break;
1610     default:
1611         CV_Error(Error::StsBadArg, "Unknown comparison operation");
1612     }
1613 }
1614
1615
1616 template<typename _Tp, typename _WTp> static void
1617 compareS_(const _Tp* src1, _WTp value, uchar* dst, size_t total, int cmpop)
1618 {
1619     size_t i;
1620     switch( cmpop )
1621     {
1622     case CMP_LT:
1623         for( i = 0; i < total; i++ )
1624             dst[i] = src1[i] < value ? 255 : 0;
1625         break;
1626     case CMP_LE:
1627         for( i = 0; i < total; i++ )
1628             dst[i] = src1[i] <= value ? 255 : 0;
1629         break;
1630     case CMP_EQ:
1631         for( i = 0; i < total; i++ )
1632             dst[i] = src1[i] == value ? 255 : 0;
1633         break;
1634     case CMP_NE:
1635         for( i = 0; i < total; i++ )
1636             dst[i] = src1[i] != value ? 255 : 0;
1637         break;
1638     case CMP_GE:
1639         for( i = 0; i < total; i++ )
1640             dst[i] = src1[i] >= value ? 255 : 0;
1641         break;
1642     case CMP_GT:
1643         for( i = 0; i < total; i++ )
1644             dst[i] = src1[i] > value ? 255 : 0;
1645         break;
1646     default:
1647         CV_Error(Error::StsBadArg, "Unknown comparison operation");
1648     }
1649 }
1650
1651
1652 void compare(const Mat& src1, const Mat& src2, Mat& dst, int cmpop)
1653 {
1654     CV_Assert( src1.type() == src2.type() && src1.channels() == 1 && src1.size == src2.size );
1655     dst.create( src1.dims, &src1.size[0], CV_8U );
1656     const Mat *arrays[]={&src1, &src2, &dst, 0};
1657     Mat planes[3];
1658
1659     NAryMatIterator it(arrays, planes);
1660     size_t total = planes[0].total();
1661     size_t i, nplanes = it.nplanes;
1662     int depth = src1.depth();
1663
1664     for( i = 0; i < nplanes; i++, ++it )
1665     {
1666         const uchar* sptr1 = planes[0].ptr();
1667         const uchar* sptr2 = planes[1].ptr();
1668         uchar* dptr = planes[2].ptr();
1669
1670         switch( depth )
1671         {
1672         case CV_8U:
1673             compare_((const uchar*)sptr1, (const uchar*)sptr2, dptr, total, cmpop);
1674             break;
1675         case CV_8S:
1676             compare_((const schar*)sptr1, (const schar*)sptr2, dptr, total, cmpop);
1677             break;
1678         case CV_16U:
1679             compare_((const ushort*)sptr1, (const ushort*)sptr2, dptr, total, cmpop);
1680             break;
1681         case CV_16S:
1682             compare_((const short*)sptr1, (const short*)sptr2, dptr, total, cmpop);
1683             break;
1684         case CV_32S:
1685             compare_((const int*)sptr1, (const int*)sptr2, dptr, total, cmpop);
1686             break;
1687         case CV_32F:
1688             compare_((const float*)sptr1, (const float*)sptr2, dptr, total, cmpop);
1689             break;
1690         case CV_64F:
1691             compare_((const double*)sptr1, (const double*)sptr2, dptr, total, cmpop);
1692             break;
1693         default:
1694             CV_Error(Error::StsUnsupportedFormat, "");
1695         }
1696     }
1697 }
1698
1699 void compare(const Mat& src, double value, Mat& dst, int cmpop)
1700 {
1701     CV_Assert( src.channels() == 1 );
1702     dst.create( src.dims, &src.size[0], CV_8U );
1703     const Mat *arrays[]={&src, &dst, 0};
1704     Mat planes[2];
1705
1706     NAryMatIterator it(arrays, planes);
1707     size_t total = planes[0].total();
1708     size_t i, nplanes = it.nplanes;
1709     int depth = src.depth();
1710     int ivalue = saturate_cast<int>(value);
1711
1712     for( i = 0; i < nplanes; i++, ++it )
1713     {
1714         const uchar* sptr = planes[0].ptr();
1715         uchar* dptr = planes[1].ptr();
1716
1717         switch( depth )
1718         {
1719         case CV_8U:
1720             compareS_((const uchar*)sptr, ivalue, dptr, total, cmpop);
1721             break;
1722         case CV_8S:
1723             compareS_((const schar*)sptr, ivalue, dptr, total, cmpop);
1724             break;
1725         case CV_16U:
1726             compareS_((const ushort*)sptr, ivalue, dptr, total, cmpop);
1727             break;
1728         case CV_16S:
1729             compareS_((const short*)sptr, ivalue, dptr, total, cmpop);
1730             break;
1731         case CV_32S:
1732             compareS_((const int*)sptr, ivalue, dptr, total, cmpop);
1733             break;
1734         case CV_32F:
1735             compareS_((const float*)sptr, value, dptr, total, cmpop);
1736             break;
1737         case CV_64F:
1738             compareS_((const double*)sptr, value, dptr, total, cmpop);
1739             break;
1740         default:
1741             CV_Error(Error::StsUnsupportedFormat, "");
1742         }
1743     }
1744 }
1745
1746
1747 template<typename _Tp> double
1748 cmpUlpsInt_(const _Tp* src1, const _Tp* src2, size_t total, int imaxdiff,
1749            size_t startidx, size_t& idx)
1750 {
1751     size_t i;
1752     int realmaxdiff = 0;
1753     for( i = 0; i < total; i++ )
1754     {
1755         int diff = std::abs(src1[i] - src2[i]);
1756         if( realmaxdiff < diff )
1757         {
1758             realmaxdiff = diff;
1759             if( diff > imaxdiff && idx == 0 )
1760                 idx = i + startidx;
1761         }
1762     }
1763     return realmaxdiff;
1764 }
1765
1766
1767 template<> double cmpUlpsInt_<int>(const int* src1, const int* src2,
1768                                           size_t total, int imaxdiff,
1769                                           size_t startidx, size_t& idx)
1770 {
1771     size_t i;
1772     double realmaxdiff = 0;
1773     for( i = 0; i < total; i++ )
1774     {
1775         double diff = fabs((double)src1[i] - (double)src2[i]);
1776         if( realmaxdiff < diff )
1777         {
1778             realmaxdiff = diff;
1779             if( diff > imaxdiff && idx == 0 )
1780                 idx = i + startidx;
1781         }
1782     }
1783     return realmaxdiff;
1784 }
1785
1786
1787 static double
1788 cmpUlpsFlt_(const int* src1, const int* src2, size_t total, int imaxdiff, size_t startidx, size_t& idx)
1789 {
1790     const int C = 0x7fffffff;
1791     int realmaxdiff = 0;
1792     size_t i;
1793     for( i = 0; i < total; i++ )
1794     {
1795         int a = src1[i], b = src2[i];
1796         if( a < 0 ) a ^= C;
1797         if( b < 0 ) b ^= C;
1798         int diff = std::abs(a - b);
1799         if( realmaxdiff < diff )
1800         {
1801             realmaxdiff = diff;
1802             if( diff > imaxdiff && idx == 0 )
1803                 idx = i + startidx;
1804         }
1805     }
1806     return realmaxdiff;
1807 }
1808
1809
1810 static double
1811 cmpUlpsFlt_(const int64* src1, const int64* src2, size_t total, int imaxdiff, size_t startidx, size_t& idx)
1812 {
1813     const int64 C = CV_BIG_INT(0x7fffffffffffffff);
1814     double realmaxdiff = 0;
1815     size_t i;
1816     for( i = 0; i < total; i++ )
1817     {
1818         int64 a = src1[i], b = src2[i];
1819         if( a < 0 ) a ^= C;
1820         if( b < 0 ) b ^= C;
1821         double diff = fabs((double)a - (double)b);
1822         if( realmaxdiff < diff )
1823         {
1824             realmaxdiff = diff;
1825             if( diff > imaxdiff && idx == 0 )
1826                 idx = i + startidx;
1827         }
1828     }
1829     return realmaxdiff;
1830 }
1831
1832 bool cmpUlps(const Mat& src1, const Mat& src2, int imaxDiff, double* _realmaxdiff, vector<int>* loc)
1833 {
1834     CV_Assert( src1.type() == src2.type() && src1.size == src2.size );
1835     const Mat *arrays[]={&src1, &src2, 0};
1836     Mat planes[2];
1837     NAryMatIterator it(arrays, planes);
1838     size_t total = planes[0].total()*planes[0].channels();
1839     size_t i, nplanes = it.nplanes;
1840     int depth = src1.depth();
1841     size_t startidx = 1, idx = 0;
1842     if(_realmaxdiff)
1843         *_realmaxdiff = 0;
1844
1845     for( i = 0; i < nplanes; i++, ++it, startidx += total )
1846     {
1847         const uchar* sptr1 = planes[0].ptr();
1848         const uchar* sptr2 = planes[1].ptr();
1849         double realmaxdiff = 0;
1850
1851         switch( depth )
1852         {
1853         case CV_8U:
1854             realmaxdiff = cmpUlpsInt_((const uchar*)sptr1, (const uchar*)sptr2, total, imaxDiff, startidx, idx);
1855             break;
1856         case CV_8S:
1857             realmaxdiff = cmpUlpsInt_((const schar*)sptr1, (const schar*)sptr2, total, imaxDiff, startidx, idx);
1858             break;
1859         case CV_16U:
1860             realmaxdiff = cmpUlpsInt_((const ushort*)sptr1, (const ushort*)sptr2, total, imaxDiff, startidx, idx);
1861             break;
1862         case CV_16S:
1863             realmaxdiff = cmpUlpsInt_((const short*)sptr1, (const short*)sptr2, total, imaxDiff, startidx, idx);
1864             break;
1865         case CV_32S:
1866             realmaxdiff = cmpUlpsInt_((const int*)sptr1, (const int*)sptr2, total, imaxDiff, startidx, idx);
1867             break;
1868         case CV_32F:
1869             realmaxdiff = cmpUlpsFlt_((const int*)sptr1, (const int*)sptr2, total, imaxDiff, startidx, idx);
1870             break;
1871         case CV_64F:
1872             realmaxdiff = cmpUlpsFlt_((const int64*)sptr1, (const int64*)sptr2, total, imaxDiff, startidx, idx);
1873             break;
1874         default:
1875             CV_Error(Error::StsUnsupportedFormat, "");
1876         }
1877
1878         if(_realmaxdiff)
1879             *_realmaxdiff = std::max(*_realmaxdiff, realmaxdiff);
1880     }
1881     if(idx > 0 && loc)
1882         setpos(src1, *loc, idx);
1883     return idx == 0;
1884 }
1885
1886
1887 template<typename _Tp> static void
1888 checkInt_(const _Tp* a, size_t total, int imin, int imax, size_t startidx, size_t& idx)
1889 {
1890     for( size_t i = 0; i < total; i++ )
1891     {
1892         int val = a[i];
1893         if( val < imin || val > imax )
1894         {
1895             idx = i + startidx;
1896             break;
1897         }
1898     }
1899 }
1900
1901
1902 template<typename _Tp> static void
1903 checkFlt_(const _Tp* a, size_t total, double fmin, double fmax, size_t startidx, size_t& idx)
1904 {
1905     for( size_t i = 0; i < total; i++ )
1906     {
1907         double val = a[i];
1908         if( cvIsNaN(val) || cvIsInf(val) || val < fmin || val > fmax )
1909         {
1910             idx = i + startidx;
1911             break;
1912         }
1913     }
1914 }
1915
1916
1917 // checks that the array does not have NaNs and/or Infs and all the elements are
1918 // within [min_val,max_val). idx is the index of the first "bad" element.
1919 int check( const Mat& a, double fmin, double fmax, vector<int>* _idx )
1920 {
1921     const Mat *arrays[]={&a, 0};
1922     Mat plane;
1923     NAryMatIterator it(arrays, &plane);
1924     size_t total = plane.total()*plane.channels();
1925     size_t i, nplanes = it.nplanes;
1926     int depth = a.depth();
1927     size_t startidx = 1, idx = 0;
1928     int imin = 0, imax = 0;
1929
1930     if( depth <= CV_32S )
1931     {
1932         imin = cvCeil(fmin);
1933         imax = cvFloor(fmax);
1934     }
1935
1936     for( i = 0; i < nplanes; i++, ++it, startidx += total )
1937     {
1938         const uchar* aptr = plane.ptr();
1939
1940         switch( depth )
1941         {
1942             case CV_8U:
1943                 checkInt_((const uchar*)aptr, total, imin, imax, startidx, idx);
1944                 break;
1945             case CV_8S:
1946                 checkInt_((const schar*)aptr, total, imin, imax, startidx, idx);
1947                 break;
1948             case CV_16U:
1949                 checkInt_((const ushort*)aptr, total, imin, imax, startidx, idx);
1950                 break;
1951             case CV_16S:
1952                 checkInt_((const short*)aptr, total, imin, imax, startidx, idx);
1953                 break;
1954             case CV_32S:
1955                 checkInt_((const int*)aptr, total, imin, imax, startidx, idx);
1956                 break;
1957             case CV_32F:
1958                 checkFlt_((const float*)aptr, total, fmin, fmax, startidx, idx);
1959                 break;
1960             case CV_64F:
1961                 checkFlt_((const double*)aptr, total, fmin, fmax, startidx, idx);
1962                 break;
1963             default:
1964                 CV_Error(Error::StsUnsupportedFormat, "");
1965         }
1966
1967         if( idx != 0 )
1968             break;
1969     }
1970
1971     if(idx != 0 && _idx)
1972         setpos(a, *_idx, idx);
1973     return idx == 0 ? 0 : -1;
1974 }
1975
1976 #define CMP_EPS_OK 0
1977 #define CMP_EPS_BIG_DIFF -1
1978 #define CMP_EPS_INVALID_TEST_DATA -2 // there is NaN or Inf value in test data
1979 #define CMP_EPS_INVALID_REF_DATA -3 // there is NaN or Inf value in reference data
1980
1981 // compares two arrays. max_diff is the maximum actual difference,
1982 // success_err_level is maximum allowed difference, idx is the index of the first
1983 // element for which difference is >success_err_level
1984 // (or index of element with the maximum difference)
1985 int cmpEps( const Mat& arr, const Mat& refarr, double* _realmaxdiff,
1986             double success_err_level, vector<int>* _idx,
1987             bool element_wise_relative_error )
1988 {
1989     CV_Assert( arr.type() == refarr.type() && arr.size == refarr.size );
1990
1991     int ilevel = refarr.depth() <= CV_32S ? cvFloor(success_err_level) : 0;
1992     int result = CMP_EPS_OK;
1993
1994     const Mat *arrays[]={&arr, &refarr, 0};
1995     Mat planes[2];
1996     NAryMatIterator it(arrays, planes);
1997     size_t total = planes[0].total()*planes[0].channels(), j = total;
1998     size_t i, nplanes = it.nplanes;
1999     int depth = arr.depth();
2000     size_t startidx = 1, idx = 0;
2001     double realmaxdiff = 0, maxval = 0;
2002
2003     if(_realmaxdiff)
2004         *_realmaxdiff = 0;
2005
2006     if( refarr.depth() >= CV_32F && !element_wise_relative_error )
2007     {
2008         maxval = cvtest::norm( refarr, NORM_INF );
2009         maxval = MAX(maxval, 1.);
2010     }
2011
2012     for( i = 0; i < nplanes; i++, ++it, startidx += total )
2013     {
2014         const uchar* sptr1 = planes[0].ptr();
2015         const uchar* sptr2 = planes[1].ptr();
2016
2017         switch( depth )
2018         {
2019         case CV_8U:
2020             realmaxdiff = cmpUlpsInt_((const uchar*)sptr1, (const uchar*)sptr2, total, ilevel, startidx, idx);
2021             break;
2022         case CV_8S:
2023             realmaxdiff = cmpUlpsInt_((const schar*)sptr1, (const schar*)sptr2, total, ilevel, startidx, idx);
2024             break;
2025         case CV_16U:
2026             realmaxdiff = cmpUlpsInt_((const ushort*)sptr1, (const ushort*)sptr2, total, ilevel, startidx, idx);
2027             break;
2028         case CV_16S:
2029             realmaxdiff = cmpUlpsInt_((const short*)sptr1, (const short*)sptr2, total, ilevel, startidx, idx);
2030             break;
2031         case CV_32S:
2032             realmaxdiff = cmpUlpsInt_((const int*)sptr1, (const int*)sptr2, total, ilevel, startidx, idx);
2033             break;
2034         case CV_32F:
2035             for( j = 0; j < total; j++ )
2036             {
2037                 double a_val = ((float*)sptr1)[j];
2038                 double b_val = ((float*)sptr2)[j];
2039                 double threshold;
2040                 if( ((int*)sptr1)[j] == ((int*)sptr2)[j] )
2041                     continue;
2042                 if( cvIsNaN(a_val) || cvIsInf(a_val) )
2043                 {
2044                     result = CMP_EPS_INVALID_TEST_DATA;
2045                     idx = startidx + j;
2046                     break;
2047                 }
2048                 if( cvIsNaN(b_val) || cvIsInf(b_val) )
2049                 {
2050                     result = CMP_EPS_INVALID_REF_DATA;
2051                     idx = startidx + j;
2052                     break;
2053                 }
2054                 a_val = fabs(a_val - b_val);
2055                 threshold = element_wise_relative_error ? fabs(b_val) + 1 : maxval;
2056                 if( a_val > threshold*success_err_level )
2057                 {
2058                     realmaxdiff = a_val/threshold;
2059                     if( idx == 0 )
2060                         idx = startidx + j;
2061                     break;
2062                 }
2063             }
2064             break;
2065         case CV_64F:
2066             for( j = 0; j < total; j++ )
2067             {
2068                 double a_val = ((double*)sptr1)[j];
2069                 double b_val = ((double*)sptr2)[j];
2070                 double threshold;
2071                 if( ((int64*)sptr1)[j] == ((int64*)sptr2)[j] )
2072                     continue;
2073                 if( cvIsNaN(a_val) || cvIsInf(a_val) )
2074                 {
2075                     result = CMP_EPS_INVALID_TEST_DATA;
2076                     idx = startidx + j;
2077                     break;
2078                 }
2079                 if( cvIsNaN(b_val) || cvIsInf(b_val) )
2080                 {
2081                     result = CMP_EPS_INVALID_REF_DATA;
2082                     idx = startidx + j;
2083                     break;
2084                 }
2085                 a_val = fabs(a_val - b_val);
2086                 threshold = element_wise_relative_error ? fabs(b_val) + 1 : maxval;
2087                 if( a_val > threshold*success_err_level )
2088                 {
2089                     realmaxdiff = a_val/threshold;
2090                     idx = startidx + j;
2091                     break;
2092                 }
2093             }
2094             break;
2095         default:
2096             assert(0);
2097             return CMP_EPS_BIG_DIFF;
2098         }
2099         if(_realmaxdiff)
2100             *_realmaxdiff = MAX(*_realmaxdiff, realmaxdiff);
2101         if( idx != 0 )
2102             break;
2103     }
2104
2105     if( result == 0 && idx != 0 )
2106         result = CMP_EPS_BIG_DIFF;
2107
2108     if( result < -1 && _realmaxdiff )
2109         *_realmaxdiff = exp(1000.);
2110     if(idx > 0 && _idx)
2111         setpos(arr, *_idx, idx);
2112
2113     return result;
2114 }
2115
2116
2117 int cmpEps2( TS* ts, const Mat& a, const Mat& b, double success_err_level,
2118              bool element_wise_relative_error, const char* desc )
2119 {
2120     char msg[100];
2121     double diff = 0;
2122     vector<int> idx;
2123     int code = cmpEps( a, b, &diff, success_err_level, &idx, element_wise_relative_error );
2124
2125     switch( code )
2126     {
2127     case CMP_EPS_BIG_DIFF:
2128         sprintf( msg, "%s: Too big difference (=%g > %g)", desc, diff, success_err_level );
2129         code = TS::FAIL_BAD_ACCURACY;
2130         break;
2131     case CMP_EPS_INVALID_TEST_DATA:
2132         sprintf( msg, "%s: Invalid output", desc );
2133         code = TS::FAIL_INVALID_OUTPUT;
2134         break;
2135     case CMP_EPS_INVALID_REF_DATA:
2136         sprintf( msg, "%s: Invalid reference output", desc );
2137         code = TS::FAIL_INVALID_OUTPUT;
2138         break;
2139     default:
2140         ;
2141     }
2142
2143     if( code < 0 )
2144     {
2145         if( a.total() == 1 )
2146         {
2147             ts->printf( TS::LOG, "%s\n", msg );
2148         }
2149         else if( a.dims == 2 && (a.rows == 1 || a.cols == 1) )
2150         {
2151             ts->printf( TS::LOG, "%s at element %d\n", msg, idx[0] + idx[1] );
2152         }
2153         else
2154         {
2155             string idxstr = vec2str(", ", &idx[0], idx.size());
2156             ts->printf( TS::LOG, "%s at (%s)\n", msg, idxstr.c_str() );
2157         }
2158     }
2159
2160     return code;
2161 }
2162
2163
2164 int cmpEps2_64f( TS* ts, const double* val, const double* refval, int len,
2165              double eps, const char* param_name )
2166 {
2167     Mat _val(1, len, CV_64F, (void*)val);
2168     Mat _refval(1, len, CV_64F, (void*)refval);
2169
2170     return cmpEps2( ts, _val, _refval, eps, true, param_name );
2171 }
2172
2173
2174 template<typename _Tp> static void
2175 GEMM_(const _Tp* a_data0, int a_step, int a_delta,
2176       const _Tp* b_data0, int b_step, int b_delta,
2177       const _Tp* c_data0, int c_step, int c_delta,
2178       _Tp* d_data, int d_step,
2179       int d_rows, int d_cols, int a_cols, int cn,
2180       double alpha, double beta)
2181 {
2182     for( int i = 0; i < d_rows; i++, d_data += d_step, c_data0 += c_step, a_data0 += a_step )
2183     {
2184         for( int j = 0; j < d_cols; j++ )
2185         {
2186             const _Tp* a_data = a_data0;
2187             const _Tp* b_data = b_data0 + j*b_delta;
2188             const _Tp* c_data = c_data0 + j*c_delta;
2189
2190             if( cn == 1 )
2191             {
2192                 double s = 0;
2193                 for( int k = 0; k < a_cols; k++ )
2194                 {
2195                     s += ((double)a_data[0])*b_data[0];
2196                     a_data += a_delta;
2197                     b_data += b_step;
2198                 }
2199                 d_data[j] = (_Tp)(s*alpha + (c_data ? c_data[0]*beta : 0));
2200             }
2201             else
2202             {
2203                 double s_re = 0, s_im = 0;
2204
2205                 for( int k = 0; k < a_cols; k++ )
2206                 {
2207                     s_re += ((double)a_data[0])*b_data[0] - ((double)a_data[1])*b_data[1];
2208                     s_im += ((double)a_data[0])*b_data[1] + ((double)a_data[1])*b_data[0];
2209                     a_data += a_delta;
2210                     b_data += b_step;
2211                 }
2212
2213                 s_re *= alpha;
2214                 s_im *= alpha;
2215
2216                 if( c_data )
2217                 {
2218                     s_re += c_data[0]*beta;
2219                     s_im += c_data[1]*beta;
2220                 }
2221
2222                 d_data[j*2] = (_Tp)s_re;
2223                 d_data[j*2+1] = (_Tp)s_im;
2224             }
2225         }
2226     }
2227 }
2228
2229
2230 void gemm( const Mat& _a, const Mat& _b, double alpha,
2231            const Mat& _c, double beta, Mat& d, int flags )
2232 {
2233     Mat a = _a, b = _b, c = _c;
2234
2235     if( a.data == d.data )
2236         a = a.clone();
2237
2238     if( b.data == d.data )
2239         b = b.clone();
2240
2241     if( !c.empty() && c.data == d.data && (flags & cv::GEMM_3_T) )
2242         c = c.clone();
2243
2244     int a_rows = a.rows, a_cols = a.cols, b_rows = b.rows, b_cols = b.cols;
2245     int cn = a.channels();
2246     int a_step = (int)a.step1(), a_delta = cn;
2247     int b_step = (int)b.step1(), b_delta = cn;
2248     int c_rows = 0, c_cols = 0, c_step = 0, c_delta = 0;
2249
2250     CV_Assert( a.type() == b.type() && a.dims == 2 && b.dims == 2 && cn <= 2 );
2251
2252     if( flags & cv::GEMM_1_T )
2253     {
2254         std::swap( a_rows, a_cols );
2255         std::swap( a_step, a_delta );
2256     }
2257
2258     if( flags & cv::GEMM_2_T )
2259     {
2260         std::swap( b_rows, b_cols );
2261         std::swap( b_step, b_delta );
2262     }
2263
2264     if( !c.empty() )
2265     {
2266         c_rows = c.rows;
2267         c_cols = c.cols;
2268         c_step = (int)c.step1();
2269         c_delta = cn;
2270
2271         if( flags & cv::GEMM_3_T )
2272         {
2273             std::swap( c_rows, c_cols );
2274             std::swap( c_step, c_delta );
2275         }
2276
2277         CV_Assert( c.dims == 2 && c.type() == a.type() && c_rows == a_rows && c_cols == b_cols );
2278     }
2279
2280     d.create(a_rows, b_cols, a.type());
2281
2282     if( a.depth() == CV_32F )
2283         GEMM_(a.ptr<float>(), a_step, a_delta, b.ptr<float>(), b_step, b_delta,
2284               !c.empty() ? c.ptr<float>() : 0, c_step, c_delta, d.ptr<float>(),
2285               (int)d.step1(), a_rows, b_cols, a_cols, cn, alpha, beta );
2286     else
2287         GEMM_(a.ptr<double>(), a_step, a_delta, b.ptr<double>(), b_step, b_delta,
2288               !c.empty() ? c.ptr<double>() : 0, c_step, c_delta, d.ptr<double>(),
2289               (int)d.step1(), a_rows, b_cols, a_cols, cn, alpha, beta );
2290 }
2291
2292
2293 template<typename _Tp> static void
2294 transform_(const _Tp* sptr, _Tp* dptr, size_t total, int scn, int dcn, const double* mat)
2295 {
2296     for( size_t i = 0; i < total; i++, sptr += scn, dptr += dcn )
2297     {
2298         for( int j = 0; j < dcn; j++ )
2299         {
2300             double s = mat[j*(scn + 1) + scn];
2301             for( int k = 0; k < scn; k++ )
2302                 s += mat[j*(scn + 1) + k]*sptr[k];
2303             dptr[j] = saturate_cast<_Tp>(s);
2304         }
2305     }
2306 }
2307
2308
2309 void transform( const Mat& src, Mat& dst, const Mat& transmat, const Mat& _shift )
2310 {
2311     double mat[20];
2312
2313     int scn = src.channels();
2314     int dcn = dst.channels();
2315     int depth = src.depth();
2316     int mattype = transmat.depth();
2317     Mat shift = _shift.reshape(1, 0);
2318     bool haveShift = !shift.empty();
2319
2320     CV_Assert( scn <= 4 && dcn <= 4 &&
2321               (mattype == CV_32F || mattype == CV_64F) &&
2322               (!haveShift || (shift.type() == mattype && (shift.rows == 1 || shift.cols == 1))) );
2323
2324     // prepare cn x (cn + 1) transform matrix
2325     if( mattype == CV_32F )
2326     {
2327         for( int i = 0; i < transmat.rows; i++ )
2328         {
2329             mat[i*(scn+1)+scn] = 0.;
2330             for( int j = 0; j < transmat.cols; j++ )
2331                 mat[i*(scn+1)+j] = transmat.at<float>(i,j);
2332             if( haveShift )
2333                 mat[i*(scn+1)+scn] = shift.at<float>(i);
2334         }
2335     }
2336     else
2337     {
2338         for( int i = 0; i < transmat.rows; i++ )
2339         {
2340             mat[i*(scn+1)+scn] = 0.;
2341             for( int j = 0; j < transmat.cols; j++ )
2342                 mat[i*(scn+1)+j] = transmat.at<double>(i,j);
2343             if( haveShift )
2344                 mat[i*(scn+1)+scn] = shift.at<double>(i);
2345         }
2346     }
2347
2348     const Mat *arrays[]={&src, &dst, 0};
2349     Mat planes[2];
2350     NAryMatIterator it(arrays, planes);
2351     size_t total = planes[0].total();
2352     size_t i, nplanes = it.nplanes;
2353
2354     for( i = 0; i < nplanes; i++, ++it )
2355     {
2356         const uchar* sptr = planes[0].ptr();
2357         uchar* dptr = planes[1].ptr();
2358
2359         switch( depth )
2360         {
2361         case CV_8U:
2362             transform_((const uchar*)sptr, (uchar*)dptr, total, scn, dcn, mat);
2363             break;
2364         case CV_8S:
2365             transform_((const schar*)sptr, (schar*)dptr, total, scn, dcn, mat);
2366             break;
2367         case CV_16U:
2368             transform_((const ushort*)sptr, (ushort*)dptr, total, scn, dcn, mat);
2369             break;
2370         case CV_16S:
2371             transform_((const short*)sptr, (short*)dptr, total, scn, dcn, mat);
2372             break;
2373         case CV_32S:
2374             transform_((const int*)sptr, (int*)dptr, total, scn, dcn, mat);
2375             break;
2376         case CV_32F:
2377             transform_((const float*)sptr, (float*)dptr, total, scn, dcn, mat);
2378             break;
2379         case CV_64F:
2380             transform_((const double*)sptr, (double*)dptr, total, scn, dcn, mat);
2381             break;
2382         default:
2383             CV_Error(Error::StsUnsupportedFormat, "");
2384         }
2385     }
2386 }
2387
2388 template<typename _Tp> static void
2389 minmax_(const _Tp* src1, const _Tp* src2, _Tp* dst, size_t total, char op)
2390 {
2391     if( op == 'M' )
2392         for( size_t i = 0; i < total; i++ )
2393             dst[i] = std::max(src1[i], src2[i]);
2394     else
2395         for( size_t i = 0; i < total; i++ )
2396             dst[i] = std::min(src1[i], src2[i]);
2397 }
2398
2399 static void minmax(const Mat& src1, const Mat& src2, Mat& dst, char op)
2400 {
2401     dst.create(src1.dims, src1.size, src1.type());
2402     CV_Assert( src1.type() == src2.type() && src1.size == src2.size );
2403     const Mat *arrays[]={&src1, &src2, &dst, 0};
2404     Mat planes[3];
2405
2406     NAryMatIterator it(arrays, planes);
2407     size_t total = planes[0].total()*planes[0].channels();
2408     size_t i, nplanes = it.nplanes, depth = src1.depth();
2409
2410     for( i = 0; i < nplanes; i++, ++it )
2411     {
2412         const uchar* sptr1 = planes[0].ptr();
2413         const uchar* sptr2 = planes[1].ptr();
2414         uchar* dptr = planes[2].ptr();
2415
2416         switch( depth )
2417         {
2418         case CV_8U:
2419             minmax_((const uchar*)sptr1, (const uchar*)sptr2, (uchar*)dptr, total, op);
2420             break;
2421         case CV_8S:
2422             minmax_((const schar*)sptr1, (const schar*)sptr2, (schar*)dptr, total, op);
2423             break;
2424         case CV_16U:
2425             minmax_((const ushort*)sptr1, (const ushort*)sptr2, (ushort*)dptr, total, op);
2426             break;
2427         case CV_16S:
2428             minmax_((const short*)sptr1, (const short*)sptr2, (short*)dptr, total, op);
2429             break;
2430         case CV_32S:
2431             minmax_((const int*)sptr1, (const int*)sptr2, (int*)dptr, total, op);
2432             break;
2433         case CV_32F:
2434             minmax_((const float*)sptr1, (const float*)sptr2, (float*)dptr, total, op);
2435             break;
2436         case CV_64F:
2437             minmax_((const double*)sptr1, (const double*)sptr2, (double*)dptr, total, op);
2438             break;
2439         default:
2440             CV_Error(Error::StsUnsupportedFormat, "");
2441         }
2442     }
2443 }
2444
2445
2446 void min(const Mat& src1, const Mat& src2, Mat& dst)
2447 {
2448     minmax( src1, src2, dst, 'm' );
2449 }
2450
2451 void max(const Mat& src1, const Mat& src2, Mat& dst)
2452 {
2453     minmax( src1, src2, dst, 'M' );
2454 }
2455
2456
2457 template<typename _Tp> static void
2458 minmax_(const _Tp* src1, _Tp val, _Tp* dst, size_t total, char op)
2459 {
2460     if( op == 'M' )
2461         for( size_t i = 0; i < total; i++ )
2462             dst[i] = std::max(src1[i], val);
2463     else
2464         for( size_t i = 0; i < total; i++ )
2465             dst[i] = std::min(src1[i], val);
2466 }
2467
2468 static void minmax(const Mat& src1, double val, Mat& dst, char op)
2469 {
2470     dst.create(src1.dims, src1.size, src1.type());
2471     const Mat *arrays[]={&src1, &dst, 0};
2472     Mat planes[2];
2473
2474     NAryMatIterator it(arrays, planes);
2475     size_t total = planes[0].total()*planes[0].channels();
2476     size_t i, nplanes = it.nplanes, depth = src1.depth();
2477     int ival = saturate_cast<int>(val);
2478
2479     for( i = 0; i < nplanes; i++, ++it )
2480     {
2481         const uchar* sptr1 = planes[0].ptr();
2482         uchar* dptr = planes[1].ptr();
2483
2484         switch( depth )
2485         {
2486         case CV_8U:
2487             minmax_((const uchar*)sptr1, saturate_cast<uchar>(ival), (uchar*)dptr, total, op);
2488             break;
2489         case CV_8S:
2490             minmax_((const schar*)sptr1, saturate_cast<schar>(ival), (schar*)dptr, total, op);
2491             break;
2492         case CV_16U:
2493             minmax_((const ushort*)sptr1, saturate_cast<ushort>(ival), (ushort*)dptr, total, op);
2494             break;
2495         case CV_16S:
2496             minmax_((const short*)sptr1, saturate_cast<short>(ival), (short*)dptr, total, op);
2497             break;
2498         case CV_32S:
2499             minmax_((const int*)sptr1, saturate_cast<int>(ival), (int*)dptr, total, op);
2500             break;
2501         case CV_32F:
2502             minmax_((const float*)sptr1, saturate_cast<float>(val), (float*)dptr, total, op);
2503             break;
2504         case CV_64F:
2505             minmax_((const double*)sptr1, saturate_cast<double>(val), (double*)dptr, total, op);
2506             break;
2507         default:
2508             CV_Error(Error::StsUnsupportedFormat, "");
2509         }
2510     }
2511 }
2512
2513
2514 void min(const Mat& src1, double val, Mat& dst)
2515 {
2516     minmax( src1, val, dst, 'm' );
2517 }
2518
2519 void max(const Mat& src1, double val, Mat& dst)
2520 {
2521     minmax( src1, val, dst, 'M' );
2522 }
2523
2524
2525 template<typename _Tp> static void
2526 muldiv_(const _Tp* src1, const _Tp* src2, _Tp* dst, size_t total, double scale, char op)
2527 {
2528     if( op == '*' )
2529         for( size_t i = 0; i < total; i++ )
2530             dst[i] = saturate_cast<_Tp>((scale*src1[i])*src2[i]);
2531     else if( src1 )
2532         for( size_t i = 0; i < total; i++ )
2533             dst[i] = src2[i] ? saturate_cast<_Tp>((scale*src1[i])/src2[i]) : 0;
2534     else
2535         for( size_t i = 0; i < total; i++ )
2536             dst[i] = src2[i] ? saturate_cast<_Tp>(scale/src2[i]) : 0;
2537 }
2538
2539 static void muldiv(const Mat& src1, const Mat& src2, Mat& dst, double scale, char op)
2540 {
2541     dst.create(src2.dims, src2.size, src2.type());
2542     CV_Assert( src1.empty() || (src1.type() == src2.type() && src1.size == src2.size) );
2543     const Mat *arrays[]={&src1, &src2, &dst, 0};
2544     Mat planes[3];
2545
2546     NAryMatIterator it(arrays, planes);
2547     size_t total = planes[1].total()*planes[1].channels();
2548     size_t i, nplanes = it.nplanes, depth = src2.depth();
2549
2550     for( i = 0; i < nplanes; i++, ++it )
2551     {
2552         const uchar* sptr1 = planes[0].ptr();
2553         const uchar* sptr2 = planes[1].ptr();
2554         uchar* dptr = planes[2].ptr();
2555
2556         switch( depth )
2557         {
2558         case CV_8U:
2559             muldiv_((const uchar*)sptr1, (const uchar*)sptr2, (uchar*)dptr, total, scale, op);
2560             break;
2561         case CV_8S:
2562             muldiv_((const schar*)sptr1, (const schar*)sptr2, (schar*)dptr, total, scale, op);
2563             break;
2564         case CV_16U:
2565             muldiv_((const ushort*)sptr1, (const ushort*)sptr2, (ushort*)dptr, total, scale, op);
2566             break;
2567         case CV_16S:
2568             muldiv_((const short*)sptr1, (const short*)sptr2, (short*)dptr, total, scale, op);
2569             break;
2570         case CV_32S:
2571             muldiv_((const int*)sptr1, (const int*)sptr2, (int*)dptr, total, scale, op);
2572             break;
2573         case CV_32F:
2574             muldiv_((const float*)sptr1, (const float*)sptr2, (float*)dptr, total, scale, op);
2575             break;
2576         case CV_64F:
2577             muldiv_((const double*)sptr1, (const double*)sptr2, (double*)dptr, total, scale, op);
2578             break;
2579         default:
2580             CV_Error(Error::StsUnsupportedFormat, "");
2581         }
2582     }
2583 }
2584
2585
2586 void multiply(const Mat& src1, const Mat& src2, Mat& dst, double scale)
2587 {
2588     muldiv( src1, src2, dst, scale, '*' );
2589 }
2590
2591 void divide(const Mat& src1, const Mat& src2, Mat& dst, double scale)
2592 {
2593     muldiv( src1, src2, dst, scale, '/' );
2594 }
2595
2596
2597 template<typename _Tp> static void
2598 mean_(const _Tp* src, const uchar* mask, size_t total, int cn, Scalar& sum, int& nz)
2599 {
2600     if( !mask )
2601     {
2602         nz += (int)total;
2603         total *= cn;
2604         for( size_t i = 0; i < total; i += cn )
2605         {
2606             for( int c = 0; c < cn; c++ )
2607                 sum[c] += src[i + c];
2608         }
2609     }
2610     else
2611     {
2612         for( size_t i = 0; i < total; i++ )
2613             if( mask[i] )
2614             {
2615                 nz++;
2616                 for( int c = 0; c < cn; c++ )
2617                     sum[c] += src[i*cn + c];
2618             }
2619     }
2620 }
2621
2622 Scalar mean(const Mat& src, const Mat& mask)
2623 {
2624     CV_Assert(mask.empty() || (mask.type() == CV_8U && mask.size == src.size));
2625     Scalar sum;
2626     int nz = 0;
2627
2628     const Mat *arrays[]={&src, &mask, 0};
2629     Mat planes[2];
2630
2631     NAryMatIterator it(arrays, planes);
2632     size_t total = planes[0].total();
2633     size_t i, nplanes = it.nplanes;
2634     int depth = src.depth(), cn = src.channels();
2635
2636     for( i = 0; i < nplanes; i++, ++it )
2637     {
2638         const uchar* sptr = planes[0].ptr();
2639         const uchar* mptr = planes[1].ptr();
2640
2641         switch( depth )
2642         {
2643         case CV_8U:
2644             mean_((const uchar*)sptr, mptr, total, cn, sum, nz);
2645             break;
2646         case CV_8S:
2647             mean_((const schar*)sptr, mptr, total, cn, sum, nz);
2648             break;
2649         case CV_16U:
2650             mean_((const ushort*)sptr, mptr, total, cn, sum, nz);
2651             break;
2652         case CV_16S:
2653             mean_((const short*)sptr, mptr, total, cn, sum, nz);
2654             break;
2655         case CV_32S:
2656             mean_((const int*)sptr, mptr, total, cn, sum, nz);
2657             break;
2658         case CV_32F:
2659             mean_((const float*)sptr, mptr, total, cn, sum, nz);
2660             break;
2661         case CV_64F:
2662             mean_((const double*)sptr, mptr, total, cn, sum, nz);
2663             break;
2664         default:
2665             CV_Error(Error::StsUnsupportedFormat, "");
2666         }
2667     }
2668
2669     return sum * (1./std::max(nz, 1));
2670 }
2671
2672
2673 void  patchZeros( Mat& mat, double level )
2674 {
2675     int j, ncols = mat.cols * mat.channels();
2676     int depth = mat.depth();
2677     CV_Assert( depth == CV_32F || depth == CV_64F );
2678
2679     for( int i = 0; i < mat.rows; i++ )
2680     {
2681         if( depth == CV_32F )
2682         {
2683             float* data = mat.ptr<float>(i);
2684             for( j = 0; j < ncols; j++ )
2685                 if( fabs(data[j]) < level )
2686                     data[j] += 1;
2687         }
2688         else
2689         {
2690             double* data = mat.ptr<double>(i);
2691             for( j = 0; j < ncols; j++ )
2692                 if( fabs(data[j]) < level )
2693                     data[j] += 1;
2694         }
2695     }
2696 }
2697
2698
2699 static void calcSobelKernel1D( int order, int _aperture_size, int size, vector<int>& kernel )
2700 {
2701     int i, j, oldval, newval;
2702     kernel.resize(size + 1);
2703
2704     if( _aperture_size < 0 )
2705     {
2706         static const int scharr[] = { 3, 10, 3, -1, 0, 1 };
2707         assert( size == 3 );
2708         for( i = 0; i < size; i++ )
2709             kernel[i] = scharr[order*3 + i];
2710         return;
2711     }
2712
2713     for( i = 1; i <= size; i++ )
2714         kernel[i] = 0;
2715     kernel[0] = 1;
2716
2717     for( i = 0; i < size - order - 1; i++ )
2718     {
2719         oldval = kernel[0];
2720         for( j = 1; j <= size; j++ )
2721         {
2722             newval = kernel[j] + kernel[j-1];
2723             kernel[j-1] = oldval;
2724             oldval = newval;
2725         }
2726     }
2727
2728     for( i = 0; i < order; i++ )
2729     {
2730         oldval = -kernel[0];
2731         for( j = 1; j <= size; j++ )
2732         {
2733             newval = kernel[j-1] - kernel[j];
2734             kernel[j-1] = oldval;
2735             oldval = newval;
2736         }
2737     }
2738 }
2739
2740
2741 Mat calcSobelKernel2D( int dx, int dy, int _aperture_size, int origin )
2742 {
2743     CV_Assert( (_aperture_size == -1 || (_aperture_size >= 1 && _aperture_size % 2 == 1)) &&
2744               dx >= 0 && dy >= 0 && dx + dy <= 3 );
2745     Size ksize = _aperture_size == -1 ? Size(3,3) : _aperture_size > 1 ?
2746         Size(_aperture_size, _aperture_size) : dx > 0 ? Size(3, 1) : Size(1, 3);
2747
2748     Mat kernel(ksize, CV_32F);
2749     vector<int> kx, ky;
2750
2751     calcSobelKernel1D( dx, _aperture_size, ksize.width, kx );
2752     calcSobelKernel1D( dy, _aperture_size, ksize.height, ky );
2753
2754     for( int i = 0; i < kernel.rows; i++ )
2755     {
2756         float ay = (float)ky[i]*(origin && (dy & 1) ? -1 : 1);
2757         for( int j = 0; j < kernel.cols; j++ )
2758             kernel.at<float>(i, j) = kx[j]*ay;
2759     }
2760
2761     return kernel;
2762 }
2763
2764
2765 Mat calcLaplaceKernel2D( int aperture_size )
2766 {
2767     int ksize = aperture_size == 1 ? 3 : aperture_size;
2768     Mat kernel(ksize, ksize, CV_32F);
2769
2770     vector<int> kx, ky;
2771
2772     calcSobelKernel1D( 2, aperture_size, ksize, kx );
2773     if( aperture_size > 1 )
2774         calcSobelKernel1D( 0, aperture_size, ksize, ky );
2775     else
2776     {
2777         ky.resize(3);
2778         ky[0] = ky[2] = 0; ky[1] = 1;
2779     }
2780
2781     for( int i = 0; i < ksize; i++ )
2782         for( int j = 0; j < ksize; j++ )
2783             kernel.at<float>(i, j) = (float)(kx[j]*ky[i] + kx[i]*ky[j]);
2784
2785     return kernel;
2786 }
2787
2788
2789 void initUndistortMap( const Mat& _a0, const Mat& _k0, Size sz, Mat& _mapx, Mat& _mapy )
2790 {
2791     _mapx.create(sz, CV_32F);
2792     _mapy.create(sz, CV_32F);
2793
2794     double a[9], k[5]={0,0,0,0,0};
2795     Mat _a(3, 3, CV_64F, a);
2796     Mat _k(_k0.rows,_k0.cols, CV_MAKETYPE(CV_64F,_k0.channels()),k);
2797     double fx, fy, cx, cy, ifx, ify, cxn, cyn;
2798
2799     _a0.convertTo(_a, CV_64F);
2800     _k0.convertTo(_k, CV_64F);
2801     fx = a[0]; fy = a[4]; cx = a[2]; cy = a[5];
2802     ifx = 1./fx; ify = 1./fy;
2803     cxn = cx;
2804     cyn = cy;
2805
2806     for( int v = 0; v < sz.height; v++ )
2807     {
2808         for( int u = 0; u < sz.width; u++ )
2809         {
2810             double x = (u - cxn)*ifx;
2811             double y = (v - cyn)*ify;
2812             double x2 = x*x, y2 = y*y;
2813             double r2 = x2 + y2;
2814             double cdist = 1 + (k[0] + (k[1] + k[4]*r2)*r2)*r2;
2815             double x1 = x*cdist + k[2]*2*x*y + k[3]*(r2 + 2*x2);
2816             double y1 = y*cdist + k[3]*2*x*y + k[2]*(r2 + 2*y2);
2817
2818             _mapy.at<float>(v, u) = (float)(y1*fy + cy);
2819             _mapx.at<float>(v, u) = (float)(x1*fx + cx);
2820         }
2821     }
2822 }
2823
2824
2825 std::ostream& operator << (std::ostream& out, const MatInfo& m)
2826 {
2827     if( !m.m || m.m->empty() )
2828         out << "<Empty>";
2829     else
2830     {
2831         static const char* depthstr[] = {"8u", "8s", "16u", "16s", "32s", "32f", "64f", "?"};
2832         out << depthstr[m.m->depth()] << "C" << m.m->channels() << " " << m.m->dims << "-dim (";
2833         for( int i = 0; i < m.m->dims; i++ )
2834             out << m.m->size[i] << (i < m.m->dims-1 ? " x " : ")");
2835     }
2836     return out;
2837 }
2838
2839
2840 static Mat getSubArray(const Mat& m, int border, vector<int>& ofs0, vector<int>& ofs)
2841 {
2842     ofs.resize(ofs0.size());
2843     if( border < 0 )
2844     {
2845         std::copy(ofs0.begin(), ofs0.end(), ofs.begin());
2846         return m;
2847     }
2848     int i, d = m.dims;
2849     CV_Assert(d == (int)ofs.size());
2850     vector<Range> r(d);
2851     for( i = 0; i < d; i++ )
2852     {
2853         r[i].start = std::max(0, ofs0[i] - border);
2854         r[i].end = std::min(ofs0[i] + 1 + border, m.size[i]);
2855         ofs[i] = std::min(ofs0[i], border);
2856     }
2857     return m(&r[0]);
2858 }
2859
2860 template<typename _Tp, typename _WTp> static void
2861 writeElems(std::ostream& out, const void* data, int nelems, int starpos)
2862 {
2863     for(int i = 0; i < nelems; i++)
2864     {
2865         if( i == starpos )
2866             out << "*";
2867         out << (_WTp)((_Tp*)data)[i];
2868         if( i == starpos )
2869             out << "*";
2870         out << (i+1 < nelems ? ", " : "");
2871     }
2872 }
2873
2874
2875 static void writeElems(std::ostream& out, const void* data, int nelems, int depth, int starpos)
2876 {
2877     if(depth == CV_8U)
2878         writeElems<uchar, int>(out, data, nelems, starpos);
2879     else if(depth == CV_8S)
2880         writeElems<schar, int>(out, data, nelems, starpos);
2881     else if(depth == CV_16U)
2882         writeElems<ushort, int>(out, data, nelems, starpos);
2883     else if(depth == CV_16S)
2884         writeElems<short, int>(out, data, nelems, starpos);
2885     else if(depth == CV_32S)
2886         writeElems<int, int>(out, data, nelems, starpos);
2887     else if(depth == CV_32F)
2888     {
2889         std::streamsize pp = out.precision();
2890         out.precision(8);
2891         writeElems<float, float>(out, data, nelems, starpos);
2892         out.precision(pp);
2893     }
2894     else if(depth == CV_64F)
2895     {
2896         std::streamsize pp = out.precision();
2897         out.precision(16);
2898         writeElems<double, double>(out, data, nelems, starpos);
2899         out.precision(pp);
2900     }
2901     else
2902         CV_Error(Error::StsUnsupportedFormat, "");
2903 }
2904
2905
2906 struct MatPart
2907 {
2908     MatPart(const Mat& _m, const vector<int>* _loc)
2909     : m(&_m), loc(_loc) {}
2910     const Mat* m;
2911     const vector<int>* loc;
2912 };
2913
2914 static std::ostream& operator << (std::ostream& out, const MatPart& m)
2915 {
2916     CV_Assert( !m.loc || ((int)m.loc->size() == m.m->dims && m.m->dims <= 2) );
2917     if( !m.loc )
2918         out << *m.m;
2919     else
2920     {
2921         int i, depth = m.m->depth(), cn = m.m->channels(), width = m.m->cols*cn;
2922         for( i = 0; i < m.m->rows; i++ )
2923         {
2924             writeElems(out, m.m->ptr(i), width, depth, i == (*m.loc)[0] ? (*m.loc)[1] : -1);
2925             out << (i < m.m->rows-1 ? ";\n" : "");
2926         }
2927     }
2928     return out;
2929 }
2930
2931 MatComparator::MatComparator(double _maxdiff, int _context)
2932     : maxdiff(_maxdiff), realmaxdiff(DBL_MAX), context(_context) {}
2933
2934 ::testing::AssertionResult
2935 MatComparator::operator()(const char* expr1, const char* expr2,
2936                           const Mat& m1, const Mat& m2)
2937 {
2938     if( m1.type() != m2.type() || m1.size != m2.size )
2939         return ::testing::AssertionFailure()
2940         << "The reference and the actual output arrays have different type or size:\n"
2941         << expr1 << " ~ " << MatInfo(m1) << "\n"
2942         << expr2 << " ~ " << MatInfo(m2) << "\n";
2943
2944     //bool ok = cvtest::cmpUlps(m1, m2, maxdiff, &realmaxdiff, &loc0);
2945     int code = cmpEps( m1, m2, &realmaxdiff, maxdiff, &loc0, true);
2946
2947     if(code >= 0)
2948         return ::testing::AssertionSuccess();
2949
2950     Mat m[] = {m1.reshape(1,0), m2.reshape(1,0)};
2951     int dims = m[0].dims;
2952     vector<int> loc;
2953     int border = dims <= 2 ? context : 0;
2954
2955     Mat m1part, m2part;
2956     if( border == 0 )
2957     {
2958         loc = loc0;
2959         m1part = Mat(1, 1, m[0].depth(), m[0].ptr(&loc[0]));
2960         m2part = Mat(1, 1, m[1].depth(), m[1].ptr(&loc[0]));
2961     }
2962     else
2963     {
2964         m1part = getSubArray(m[0], border, loc0, loc);
2965         m2part = getSubArray(m[1], border, loc0, loc);
2966     }
2967
2968     return ::testing::AssertionFailure()
2969     << "too big relative difference (" << realmaxdiff << " > "
2970     << maxdiff << ") between "
2971     << MatInfo(m1) << " '" << expr1 << "' and '" << expr2 << "' at " << Mat(loc0).t() << ".\n"
2972     << "- " << expr1 << ":\n" << MatPart(m1part, border > 0 ? &loc : 0) << ".\n"
2973     << "- " << expr2 << ":\n" << MatPart(m2part, border > 0 ? &loc : 0) << ".\n";
2974 }
2975
2976 void printVersionInfo(bool useStdOut)
2977 {
2978     // Tell CTest not to discard any output
2979     if(useStdOut) std::cout << "CTEST_FULL_OUTPUT" << std::endl;
2980
2981     ::testing::Test::RecordProperty("cv_version", CV_VERSION);
2982     if(useStdOut) std::cout << "OpenCV version: " << CV_VERSION << std::endl;
2983
2984     std::string buildInfo( cv::getBuildInformation() );
2985
2986     size_t pos1 = buildInfo.find("Version control");
2987     size_t pos2 = buildInfo.find('\n', pos1);
2988     if(pos1 != std::string::npos && pos2 != std::string::npos)
2989     {
2990         size_t value_start = buildInfo.rfind(' ', pos2) + 1;
2991         std::string ver( buildInfo.substr(value_start, pos2 - value_start) );
2992         ::testing::Test::RecordProperty("cv_vcs_version", ver);
2993         if (useStdOut) std::cout << "OpenCV VCS version: " << ver << std::endl;
2994     }
2995
2996     pos1 = buildInfo.find("inner version");
2997     pos2 = buildInfo.find('\n', pos1);
2998     if(pos1 != std::string::npos && pos2 != std::string::npos)
2999     {
3000         size_t value_start = buildInfo.rfind(' ', pos2) + 1;
3001         std::string ver( buildInfo.substr(value_start, pos2 - value_start) );
3002         ::testing::Test::RecordProperty("cv_inner_vcs_version", ver);
3003         if(useStdOut) std::cout << "Inner VCS version: " << ver << std::endl;
3004     }
3005
3006     const char * build_type =
3007 #ifdef _DEBUG
3008         "debug";
3009 #else
3010         "release";
3011 #endif
3012
3013     ::testing::Test::RecordProperty("cv_build_type", build_type);
3014     if (useStdOut) std::cout << "Build type: " << build_type << std::endl;
3015
3016     const char* parallel_framework = currentParallelFramework();
3017
3018     if (parallel_framework) {
3019         ::testing::Test::RecordProperty("cv_parallel_framework", parallel_framework);
3020         if (useStdOut) std::cout << "Parallel framework: " << parallel_framework << std::endl;
3021     }
3022
3023     std::string cpu_features;
3024
3025 #if CV_POPCNT
3026     if (checkHardwareSupport(CV_CPU_POPCNT)) cpu_features += " popcnt";
3027 #endif
3028 #if CV_MMX
3029     if (checkHardwareSupport(CV_CPU_MMX)) cpu_features += " mmx";
3030 #endif
3031 #if CV_SSE
3032     if (checkHardwareSupport(CV_CPU_SSE)) cpu_features += " sse";
3033 #endif
3034 #if CV_SSE2
3035     if (checkHardwareSupport(CV_CPU_SSE2)) cpu_features += " sse2";
3036 #endif
3037 #if CV_SSE3
3038     if (checkHardwareSupport(CV_CPU_SSE3)) cpu_features += " sse3";
3039 #endif
3040 #if CV_SSSE3
3041     if (checkHardwareSupport(CV_CPU_SSSE3)) cpu_features += " ssse3";
3042 #endif
3043 #if CV_SSE4_1
3044     if (checkHardwareSupport(CV_CPU_SSE4_1)) cpu_features += " sse4.1";
3045 #endif
3046 #if CV_SSE4_2
3047     if (checkHardwareSupport(CV_CPU_SSE4_2)) cpu_features += " sse4.2";
3048 #endif
3049 #if CV_AVX
3050     if (checkHardwareSupport(CV_CPU_AVX)) cpu_features += " avx";
3051 #endif
3052 #if CV_AVX2
3053     if (checkHardwareSupport(CV_CPU_AVX2)) cpu_features += " avx2";
3054 #endif
3055 #if CV_FMA3
3056     if (checkHardwareSupport(CV_CPU_FMA3)) cpu_features += " fma3";
3057 #endif
3058 #if CV_AVX_512F
3059     if (checkHardwareSupport(CV_CPU_AVX_512F)) cpu_features += " avx-512f";
3060 #endif
3061 #if CV_AVX_512BW
3062     if (checkHardwareSupport(CV_CPU_AVX_512BW)) cpu_features += " avx-512bw";
3063 #endif
3064 #if CV_AVX_512CD
3065     if (checkHardwareSupport(CV_CPU_AVX_512CD)) cpu_features += " avx-512cd";
3066 #endif
3067 #if CV_AVX_512DQ
3068     if (checkHardwareSupport(CV_CPU_AVX_512DQ)) cpu_features += " avx-512dq";
3069 #endif
3070 #if CV_AVX_512ER
3071     if (checkHardwareSupport(CV_CPU_AVX_512ER)) cpu_features += " avx-512er";
3072 #endif
3073 #if CV_AVX_512IFMA512
3074     if (checkHardwareSupport(CV_CPU_AVX_512IFMA512)) cpu_features += " avx-512ifma512";
3075 #endif
3076 #if CV_AVX_512PF
3077     if (checkHardwareSupport(CV_CPU_AVX_512PF)) cpu_features += " avx-512pf";
3078 #endif
3079 #if CV_AVX_512VBMI
3080     if (checkHardwareSupport(CV_CPU_AVX_512VBMI)) cpu_features += " avx-512vbmi";
3081 #endif
3082 #if CV_AVX_512VL
3083     if (checkHardwareSupport(CV_CPU_AVX_512VL)) cpu_features += " avx-512vl";
3084 #endif
3085 #if CV_NEON
3086     if (checkHardwareSupport(CV_CPU_NEON)) cpu_features += " neon";
3087 #endif
3088 #if CV_FP16
3089     if (checkHardwareSupport(CV_CPU_FP16)) cpu_features += " fp16";
3090 #endif
3091 #if CV_VSX
3092     if (checkHardwareSupport(CV_CPU_VSX)) cpu_features += " VSX";
3093 #endif
3094
3095     cpu_features.erase(0, 1); // erase initial space
3096
3097     ::testing::Test::RecordProperty("cv_cpu_features", cpu_features);
3098     if (useStdOut) std::cout << "CPU features: " << cpu_features << std::endl;
3099
3100 #ifdef HAVE_IPP
3101     const char * ipp_optimization = cv::ipp::useIPP()? "enabled" : "disabled";
3102     ::testing::Test::RecordProperty("cv_ipp_optimization", ipp_optimization);
3103     if (useStdOut) std::cout << "Intel(R) IPP optimization: " << ipp_optimization << std::endl;
3104
3105     cv::String ippVer = cv::ipp::getIppVersion();
3106     ::testing::Test::RecordProperty("cv_ipp_version", ippVer);
3107     if(useStdOut) std::cout << "Intel(R) IPP version: " << ippVer.c_str() << std::endl;
3108 #endif
3109 }
3110
3111
3112
3113 void threshold( const Mat& _src, Mat& _dst,
3114                             double thresh, double maxval, int thresh_type )
3115 {
3116     int i, j;
3117     int depth = _src.depth(), cn = _src.channels();
3118     int width_n = _src.cols*cn, height = _src.rows;
3119     int ithresh = cvFloor(thresh);
3120     int imaxval, ithresh2;
3121
3122     if( depth == CV_8U )
3123     {
3124         ithresh2 = saturate_cast<uchar>(ithresh);
3125         imaxval = saturate_cast<uchar>(maxval);
3126     }
3127     else if( depth == CV_16S )
3128     {
3129         ithresh2 = saturate_cast<short>(ithresh);
3130         imaxval = saturate_cast<short>(maxval);
3131     }
3132     else
3133     {
3134         ithresh2 = cvRound(ithresh);
3135         imaxval = cvRound(maxval);
3136     }
3137
3138     assert( depth == CV_8U || depth == CV_16S || depth == CV_32F );
3139
3140     switch( thresh_type )
3141     {
3142     case CV_THRESH_BINARY:
3143         for( i = 0; i < height; i++ )
3144         {
3145             if( depth == CV_8U )
3146             {
3147                 const uchar* src = _src.ptr<uchar>(i);
3148                 uchar* dst = _dst.ptr<uchar>(i);
3149                 for( j = 0; j < width_n; j++ )
3150                     dst[j] = (uchar)(src[j] > ithresh ? imaxval : 0);
3151             }
3152             else if( depth == CV_16S )
3153             {
3154                 const short* src = _src.ptr<short>(i);
3155                 short* dst = _dst.ptr<short>(i);
3156                 for( j = 0; j < width_n; j++ )
3157                     dst[j] = (short)(src[j] > ithresh ? imaxval : 0);
3158             }
3159             else
3160             {
3161                 const float* src = _src.ptr<float>(i);
3162                 float* dst = _dst.ptr<float>(i);
3163                 for( j = 0; j < width_n; j++ )
3164                     dst[j] = (float)((double)src[j] > thresh ? maxval : 0.f);
3165             }
3166         }
3167         break;
3168     case CV_THRESH_BINARY_INV:
3169         for( i = 0; i < height; i++ )
3170         {
3171             if( depth == CV_8U )
3172             {
3173                 const uchar* src = _src.ptr<uchar>(i);
3174                 uchar* dst = _dst.ptr<uchar>(i);
3175                 for( j = 0; j < width_n; j++ )
3176                     dst[j] = (uchar)(src[j] > ithresh ? 0 : imaxval);
3177             }
3178             else if( depth == CV_16S )
3179             {
3180                 const short* src = _src.ptr<short>(i);
3181                 short* dst = _dst.ptr<short>(i);
3182                 for( j = 0; j < width_n; j++ )
3183                     dst[j] = (short)(src[j] > ithresh ? 0 : imaxval);
3184             }
3185             else
3186             {
3187                 const float* src = _src.ptr<float>(i);
3188                 float* dst = _dst.ptr<float>(i);
3189                 for( j = 0; j < width_n; j++ )
3190                     dst[j] = (float)((double)src[j] > thresh ? 0.f : maxval);
3191             }
3192         }
3193         break;
3194     case CV_THRESH_TRUNC:
3195         for( i = 0; i < height; i++ )
3196         {
3197             if( depth == CV_8U )
3198             {
3199                 const uchar* src = _src.ptr<uchar>(i);
3200                 uchar* dst = _dst.ptr<uchar>(i);
3201                 for( j = 0; j < width_n; j++ )
3202                 {
3203                     int s = src[j];
3204                     dst[j] = (uchar)(s > ithresh ? ithresh2 : s);
3205                 }
3206             }
3207             else if( depth == CV_16S )
3208             {
3209                 const short* src = _src.ptr<short>(i);
3210                 short* dst = _dst.ptr<short>(i);
3211                 for( j = 0; j < width_n; j++ )
3212                 {
3213                     int s = src[j];
3214                     dst[j] = (short)(s > ithresh ? ithresh2 : s);
3215                 }
3216             }
3217             else
3218             {
3219                 const float* src = _src.ptr<float>(i);
3220                 float* dst = _dst.ptr<float>(i);
3221                 for( j = 0; j < width_n; j++ )
3222                 {
3223                     double s = src[j];
3224                     dst[j] = (float)(s > thresh ? thresh : s);
3225                 }
3226             }
3227         }
3228         break;
3229     case CV_THRESH_TOZERO:
3230         for( i = 0; i < height; i++ )
3231         {
3232             if( depth == CV_8U )
3233             {
3234                 const uchar* src = _src.ptr<uchar>(i);
3235                 uchar* dst = _dst.ptr<uchar>(i);
3236                 for( j = 0; j < width_n; j++ )
3237                 {
3238                     int s = src[j];
3239                     dst[j] = (uchar)(s > ithresh ? s : 0);
3240                 }
3241             }
3242             else if( depth == CV_16S )
3243             {
3244                 const short* src = _src.ptr<short>(i);
3245                 short* dst = _dst.ptr<short>(i);
3246                 for( j = 0; j < width_n; j++ )
3247                 {
3248                     int s = src[j];
3249                     dst[j] = (short)(s > ithresh ? s : 0);
3250                 }
3251             }
3252             else
3253             {
3254                 const float* src = _src.ptr<float>(i);
3255                 float* dst = _dst.ptr<float>(i);
3256                 for( j = 0; j < width_n; j++ )
3257                 {
3258                     float s = src[j];
3259                     dst[j] = s > thresh ? s : 0.f;
3260                 }
3261             }
3262         }
3263         break;
3264     case CV_THRESH_TOZERO_INV:
3265         for( i = 0; i < height; i++ )
3266         {
3267             if( depth == CV_8U )
3268             {
3269                 const uchar* src = _src.ptr<uchar>(i);
3270                 uchar* dst = _dst.ptr<uchar>(i);
3271                 for( j = 0; j < width_n; j++ )
3272                 {
3273                     int s = src[j];
3274                     dst[j] = (uchar)(s > ithresh ? 0 : s);
3275                 }
3276             }
3277             else if( depth == CV_16S )
3278             {
3279                 const short* src = _src.ptr<short>(i);
3280                 short* dst = _dst.ptr<short>(i);
3281                 for( j = 0; j < width_n; j++ )
3282                 {
3283                     int s = src[j];
3284                     dst[j] = (short)(s > ithresh ? 0 : s);
3285                 }
3286             }
3287             else
3288             {
3289                 const float* src = _src.ptr<float>(i);
3290                 float* dst = _dst.ptr<float>(i);
3291                 for( j = 0; j < width_n; j++ )
3292                 {
3293                     float s = src[j];
3294                     dst[j] = s > thresh ? 0.f : s;
3295                 }
3296             }
3297         }
3298         break;
3299     default:
3300         assert(0);
3301     }
3302 }
3303
3304
3305 static void
3306 _minMaxIdx( const float* src, const uchar* mask, double* _minVal, double* _maxVal,
3307             size_t* _minIdx, size_t* _maxIdx, int len, size_t startIdx )
3308 {
3309     double minVal = FLT_MAX, maxVal = -FLT_MAX;
3310     size_t minIdx = 0, maxIdx = 0;
3311
3312     if( !mask )
3313     {
3314         for( int i = 0; i < len; i++ )
3315         {
3316             float val = src[i];
3317             if( val < minVal )
3318             {
3319                 minVal = val;
3320                 minIdx = startIdx + i;
3321             }
3322             if( val > maxVal )
3323             {
3324                 maxVal = val;
3325                 maxIdx = startIdx + i;
3326             }
3327         }
3328     }
3329     else
3330     {
3331         for( int i = 0; i < len; i++ )
3332         {
3333             float val = src[i];
3334             if( mask[i] && val < minVal )
3335             {
3336                 minVal = val;
3337                 minIdx = startIdx + i;
3338             }
3339             if( mask[i] && val > maxVal )
3340             {
3341                 maxVal = val;
3342                 maxIdx = startIdx + i;
3343             }
3344         }
3345     }
3346
3347     if (_minIdx)
3348         *_minIdx = minIdx;
3349     if (_maxIdx)
3350         *_maxIdx = maxIdx;
3351     if (_minVal)
3352         *_minVal = minVal;
3353     if (_maxVal)
3354         *_maxVal = maxVal;
3355 }
3356
3357
3358 void minMaxIdx( InputArray _img, double* minVal, double* maxVal,
3359                     Point* minLoc, Point* maxLoc, InputArray _mask )
3360 {
3361     Mat img = _img.getMat();
3362     Mat mask = _mask.getMat();
3363     CV_Assert(img.dims <= 2);
3364
3365     _minMaxIdx((const float*)img.data, mask.data, minVal, maxVal, (size_t*)minLoc, (size_t*)maxLoc, (int)img.total(),1);
3366     if( minLoc )
3367         std::swap(minLoc->x, minLoc->y);
3368     if( maxLoc )
3369         std::swap(maxLoc->x, maxLoc->y);
3370 }
3371
3372 }