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