modules/core/test/test_math.cpp

   1 //////////////////////////////////////////////////////////////////////////////////////////
   2 /////////////////// tests for matrix operations and math functions ///////////////////////
   3 //////////////////////////////////////////////////////////////////////////////////////////
   4
   5 #include "test_precomp.hpp"
   6 #include <float.h>
   7 #include <math.h>
   8
   9 using namespace cv;
  10 using namespace std;
  11
  12 /// !!! NOTE !!! These tests happily avoid overflow cases & out-of-range arguments
  13 /// so that output arrays contain neigher Inf's nor Nan's.
  14 /// Handling such cases would require special modification of check function
  15 /// (validate_test_results) => TBD.
  16 /// Also, need some logarithmic-scale generation of input data. Right now it is done (in some tests)
  17 /// by generating min/max boundaries for random data in logarimithic scale, but
  18 /// within the same test case all the input array elements are of the same order.
  19
  20 class Core_MathTest : public cvtest::ArrayTest
  21 {
  22 public:
  23     typedef cvtest::ArrayTest Base;
  24     Core_MathTest();
  25 protected:
  26     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes,
  27                                         vector<vector<int> >& types);
  28     double get_success_error_level( int /*test_case_idx*/, int i, int j );
  29     bool test_nd;
  30 };
  31
  32
  33 Core_MathTest::Core_MathTest()
  34 {
  35     optional_mask = false;
  36
  37     test_array[INPUT].push_back(NULL);
  38     test_array[OUTPUT].push_back(NULL);
  39     test_array[REF_OUTPUT].push_back(NULL);
  40
  41     test_nd = false;
  42 }
  43
  44
  45 double Core_MathTest::get_success_error_level( int /*test_case_idx*/, int i, int j )
  46 {
  47     return test_mat[i][j].depth() == CV_32F ? FLT_EPSILON*128 : DBL_EPSILON*1024;
  48 }
  49
  50
  51 void Core_MathTest::get_test_array_types_and_sizes( int test_case_idx,
  52                                                      vector<vector<Size> >& sizes,
  53                                                      vector<vector<int> >& types)
  54 {
  55     RNG& rng = ts->get_rng();
  56     int depth = cvtest::randInt(rng)%2 + CV_32F;
  57     int cn = cvtest::randInt(rng) % 4 + 1, type = CV_MAKETYPE(depth, cn);
  58     size_t i, j;
  59     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
  60
  61     for( i = 0; i < test_array.size(); i++ )
  62     {
  63         size_t count = test_array[i].size();
  64         for( j = 0; j < count; j++ )
  65             types[i][j] = type;
  66     }
  67     test_nd = cvtest::randInt(rng)%3 == 0;
  68 }
  69
  70
  71 ////////// pow /////////////
  72
  73 class Core_PowTest : public Core_MathTest
  74 {
  75 public:
  76     typedef Core_MathTest Base;
  77     Core_PowTest();
  78 protected:
  79     void get_test_array_types_and_sizes( int test_case_idx,
  80                                         vector<vector<Size> >& sizes,
  81                                         vector<vector<int> >& types );
  82     void get_minmax_bounds( int i, int j, int type, Scalar& low, Scalar& high );
  83     void run_func();
  84     void prepare_to_validation( int test_case_idx );
  85     double get_success_error_level( int test_case_idx, int i, int j );
  86     double power;
  87 };
  88
  89
  90 Core_PowTest::Core_PowTest()
  91 {
  92     power = 0;
  93 }
  94
  95
  96 void Core_PowTest::get_test_array_types_and_sizes( int test_case_idx,
  97                                                     vector<vector<Size> >& sizes,
  98                                                     vector<vector<int> >& types )
  99 {
 100     RNG& rng = ts->get_rng();
 101     int depth = cvtest::randInt(rng) % (CV_64F+1);
 102     int cn = cvtest::randInt(rng) % 4 + 1;
 103     size_t i, j;
 104     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
 105     depth += depth == CV_8S;
 106
 107     if( depth < CV_32F || cvtest::randInt(rng)%8 == 0 )
 108         // integer power
 109         power = (int)(cvtest::randInt(rng)%21 - 10);
 110     else
 111     {
 112         i = cvtest::randInt(rng)%17;
 113         power = i == 16 ? 1./3 : i == 15 ? 0.5 : i == 14 ? -0.5 : cvtest::randReal(rng)*10 - 5;
 114     }
 115
 116     for( i = 0; i < test_array.size(); i++ )
 117     {
 118         size_t count = test_array[i].size();
 119         int type = CV_MAKETYPE(depth, cn);
 120         for( j = 0; j < count; j++ )
 121             types[i][j] = type;
 122     }
 123     test_nd = cvtest::randInt(rng)%3 == 0;
 124 }
 125
 126
 127 double Core_PowTest::get_success_error_level( int test_case_idx, int i, int j )
 128 {
 129     int depth = test_mat[i][j].depth();
 130     if( depth < CV_32F )
 131         return power == cvRound(power) && power >= 0 ? 0 : 1;
 132     else
 133         return Base::get_success_error_level( test_case_idx, i, j );
 134 }
 135
 136
 137 void Core_PowTest::get_minmax_bounds( int /*i*/, int /*j*/, int type, Scalar& low, Scalar& high )
 138 {
 139     double l, u = cvtest::randInt(ts->get_rng())%1000 + 1;
 140     if( power > 0 )
 141     {
 142         double mval = cvtest::getMaxVal(type);
 143         double u1 = pow(mval,1./power)*2;
 144         u = MIN(u,u1);
 145     }
 146
 147     l = power == cvRound(power) ? -u : FLT_EPSILON;
 148     low = Scalar::all(l);
 149     high = Scalar::all(u);
 150 }
 151
 152
 153 void Core_PowTest::run_func()
 154 {
 155     if(!test_nd)
 156     {
 157         if( fabs(power-1./3) <= DBL_EPSILON && test_mat[INPUT][0].depth() == CV_32F )
 158         {
 159             Mat a = test_mat[INPUT][0], b = test_mat[OUTPUT][0];
 160
 161             a = a.reshape(1);
 162             b = b.reshape(1);
 163             for( int i = 0; i < a.rows; i++ )
 164             {
 165                 b.at<float>(i,0) = (float)fabs(cvCbrt(a.at<float>(i,0)));
 166                 for( int j = 1; j < a.cols; j++ )
 167                     b.at<float>(i,j) = (float)fabs(cv::cubeRoot(a.at<float>(i,j)));
 168             }
 169         }
 170         else
 171             cvPow( test_array[INPUT][0], test_array[OUTPUT][0], power );
 172     }
 173     else
 174     {
 175         Mat& a = test_mat[INPUT][0];
 176         Mat& b = test_mat[OUTPUT][0];
 177         if(power == 0.5)
 178             cv::sqrt(a, b);
 179         else
 180             cv::pow(a, power, b);
 181     }
 182 }
 183
 184
 185 inline static int ipow( int a, int power )
 186 {
 187     int b = 1;
 188     while( power > 0 )
 189     {
 190         if( power&1 )
 191             b *= a, power--;
 192         else
 193             a *= a, power >>= 1;
 194     }
 195     return b;
 196 }
 197
 198
 199 inline static double ipow( double a, int power )
 200 {
 201     double b = 1.;
 202     while( power > 0 )
 203     {
 204         if( power&1 )
 205             b *= a, power--;
 206         else
 207             a *= a, power >>= 1;
 208     }
 209     return b;
 210 }
 211
 212
 213 void Core_PowTest::prepare_to_validation( int /*test_case_idx*/ )
 214 {
 215     const Mat& a = test_mat[INPUT][0];
 216     Mat& b = test_mat[REF_OUTPUT][0];
 217
 218     int depth = a.depth();
 219     int ncols = a.cols*a.channels();
 220     int ipower = cvRound(power), apower = abs(ipower);
 221     int i, j;
 222
 223     for( i = 0; i < a.rows; i++ )
 224     {
 225         const uchar* a_data = a.ptr(i);
 226         uchar* b_data = b.ptr(i);
 227
 228         switch( depth )
 229         {
 230             case CV_8U:
 231                 if( ipower < 0 )
 232                     for( j = 0; j < ncols; j++ )
 233                     {
 234                         int val = ((uchar*)a_data)[j];
 235                         ((uchar*)b_data)[j] = (uchar)(val <= 1 ? val :
 236                                                       val == 2 && ipower == -1 ? 1 : 0);
 237                     }
 238                 else
 239                     for( j = 0; j < ncols; j++ )
 240                     {
 241                         int val = ((uchar*)a_data)[j];
 242                         val = ipow( val, ipower );
 243                         ((uchar*)b_data)[j] = saturate_cast<uchar>(val);
 244                     }
 245                 break;
 246             case CV_8S:
 247                 if( ipower < 0 )
 248                     for( j = 0; j < ncols; j++ )
 249                     {
 250                         int val = ((char*)a_data)[j];
 251                         ((char*)b_data)[j] = (char)((val&~1)==0 ? val :
 252                                                     val ==-1 ? 1-2*(ipower&1) :
 253                                                     val == 2 && ipower == -1 ? 1 : 0);
 254                     }
 255                 else
 256                     for( j = 0; j < ncols; j++ )
 257                     {
 258                         int val = ((char*)a_data)[j];
 259                         val = ipow( val, ipower );
 260                         ((char*)b_data)[j] = saturate_cast<schar>(val);
 261                     }
 262                 break;
 263             case CV_16U:
 264                 if( ipower < 0 )
 265                     for( j = 0; j < ncols; j++ )
 266                     {
 267                         int val = ((ushort*)a_data)[j];
 268                         ((ushort*)b_data)[j] = (ushort)((val&~1)==0 ? val :
 269                                                         val ==-1 ? 1-2*(ipower&1) :
 270                                                         val == 2 && ipower == -1 ? 1 : 0);
 271                     }
 272                 else
 273                     for( j = 0; j < ncols; j++ )
 274                     {
 275                         int val = ((ushort*)a_data)[j];
 276                         val = ipow( val, ipower );
 277                         ((ushort*)b_data)[j] = saturate_cast<ushort>(val);
 278                     }
 279                 break;
 280             case CV_16S:
 281                 if( ipower < 0 )
 282                     for( j = 0; j < ncols; j++ )
 283                     {
 284                         int val = ((short*)a_data)[j];
 285                         ((short*)b_data)[j] = (short)((val&~1)==0 ? val :
 286                                                       val ==-1 ? 1-2*(ipower&1) :
 287                                                       val == 2 && ipower == -1 ? 1 : 0);
 288                     }
 289                 else
 290                     for( j = 0; j < ncols; j++ )
 291                     {
 292                         int val = ((short*)a_data)[j];
 293                         val = ipow( val, ipower );
 294                         ((short*)b_data)[j] = saturate_cast<short>(val);
 295                     }
 296                 break;
 297             case CV_32S:
 298                 if( ipower < 0 )
 299                     for( j = 0; j < ncols; j++ )
 300                     {
 301                         int val = ((int*)a_data)[j];
 302                         ((int*)b_data)[j] = (val&~1)==0 ? val :
 303                         val ==-1 ? 1-2*(ipower&1) :
 304                         val == 2 && ipower == -1 ? 1 : 0;
 305                     }
 306                 else
 307                     for( j = 0; j < ncols; j++ )
 308                     {
 309                         int val = ((int*)a_data)[j];
 310                         val = ipow( val, ipower );
 311                         ((int*)b_data)[j] = val;
 312                     }
 313                 break;
 314             case CV_32F:
 315                 if( power != ipower )
 316                     for( j = 0; j < ncols; j++ )
 317                     {
 318                         double val = ((float*)a_data)[j];
 319                         val = pow( fabs(val), power );
 320                         ((float*)b_data)[j] = (float)val;
 321                     }
 322                 else
 323                     for( j = 0; j < ncols; j++ )
 324                     {
 325                         double val = ((float*)a_data)[j];
 326                         if( ipower < 0 )
 327                             val = 1./val;
 328                         val = ipow( val, apower );
 329                         ((float*)b_data)[j] = (float)val;
 330                     }
 331                 break;
 332             case CV_64F:
 333                 if( power != ipower )
 334                     for( j = 0; j < ncols; j++ )
 335                     {
 336                         double val = ((double*)a_data)[j];
 337                         val = pow( fabs(val), power );
 338                         ((double*)b_data)[j] = (double)val;
 339                     }
 340                 else
 341                     for( j = 0; j < ncols; j++ )
 342                     {
 343                         double val = ((double*)a_data)[j];
 344                         if( ipower < 0 )
 345                             val = 1./val;
 346                         val = ipow( val, apower );
 347                         ((double*)b_data)[j] = (double)val;
 348                     }
 349                 break;
 350         }
 351     }
 352 }
 353
 354
 355
 356 ///////////////////////////////////////// matrix tests ////////////////////////////////////////////
 357
 358 class Core_MatrixTest : public cvtest::ArrayTest
 359 {
 360 public:
 361     typedef cvtest::ArrayTest Base;
 362     Core_MatrixTest( int in_count, int out_count,
 363                        bool allow_int, bool scalar_output, int max_cn );
 364 protected:
 365     void get_test_array_types_and_sizes( int test_case_idx,
 366                                         vector<vector<Size> >& sizes,
 367                                         vector<vector<int> >& types );
 368     double get_success_error_level( int test_case_idx, int i, int j );
 369     bool allow_int;
 370     bool scalar_output;
 371     int max_cn;
 372 };
 373
 374
 375 Core_MatrixTest::Core_MatrixTest( int in_count, int out_count,
 376                                       bool _allow_int, bool _scalar_output, int _max_cn )
 377 : allow_int(_allow_int), scalar_output(_scalar_output), max_cn(_max_cn)
 378 {
 379     int i;
 380     for( i = 0; i < in_count; i++ )
 381         test_array[INPUT].push_back(NULL);
 382
 383     for( i = 0; i < out_count; i++ )
 384     {
 385         test_array[OUTPUT].push_back(NULL);
 386         test_array[REF_OUTPUT].push_back(NULL);
 387     }
 388
 389     element_wise_relative_error = false;
 390 }
 391
 392
 393 void Core_MatrixTest::get_test_array_types_and_sizes( int test_case_idx,
 394                                                        vector<vector<Size> >& sizes,
 395                                                        vector<vector<int> >& types )
 396 {
 397     RNG& rng = ts->get_rng();
 398     int depth = cvtest::randInt(rng) % (allow_int ? CV_64F+1 : 2);
 399     int cn = cvtest::randInt(rng) % max_cn + 1;
 400     size_t i, j;
 401
 402     if( allow_int )
 403         depth += depth == CV_8S;
 404     else
 405         depth += CV_32F;
 406
 407     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
 408
 409     for( i = 0; i < test_array.size(); i++ )
 410     {
 411         size_t count = test_array[i].size();
 412         int flag = (i == OUTPUT || i == REF_OUTPUT) && scalar_output;
 413         int type = !flag ? CV_MAKETYPE(depth, cn) : CV_64FC1;
 414
 415         for( j = 0; j < count; j++ )
 416         {
 417             types[i][j] = type;
 418             if( flag )
 419                 sizes[i][j] = Size( 4, 1 );
 420         }
 421     }
 422 }
 423
 424
 425 double Core_MatrixTest::get_success_error_level( int test_case_idx, int i, int j )
 426 {
 427     int input_depth = test_mat[INPUT][0].depth();
 428     double input_precision = input_depth < CV_32F ? 0 : input_depth == CV_32F ? 5e-5 : 5e-10;
 429     double output_precision = Base::get_success_error_level( test_case_idx, i, j );
 430     return MAX(input_precision, output_precision);
 431 }
 432
 433
 434 ///////////////// Trace /////////////////////
 435
 436 class Core_TraceTest : public Core_MatrixTest
 437 {
 438 public:
 439     Core_TraceTest();
 440 protected:
 441     void run_func();
 442     void prepare_to_validation( int test_case_idx );
 443 };
 444
 445
 446 Core_TraceTest::Core_TraceTest() : Core_MatrixTest( 1, 1, true, true, 4 )
 447 {
 448 }
 449
 450
 451 void Core_TraceTest::run_func()
 452 {
 453     test_mat[OUTPUT][0].at<Scalar>(0,0) = cvTrace(test_array[INPUT][0]);
 454 }
 455
 456
 457 void Core_TraceTest::prepare_to_validation( int )
 458 {
 459     Mat& mat = test_mat[INPUT][0];
 460     int count = MIN( mat.rows, mat.cols );
 461     Mat diag(count, 1, mat.type(), mat.data, mat.step + mat.elemSize());
 462     Scalar r = cvtest::mean(diag);
 463     r *= (double)count;
 464
 465     test_mat[REF_OUTPUT][0].at<Scalar>(0,0) = r;
 466 }
 467
 468
 469 ///////// dotproduct //////////
 470
 471 class Core_DotProductTest : public Core_MatrixTest
 472 {
 473 public:
 474     Core_DotProductTest();
 475 protected:
 476     void run_func();
 477     void prepare_to_validation( int test_case_idx );
 478 };
 479
 480
 481 Core_DotProductTest::Core_DotProductTest() : Core_MatrixTest( 2, 1, true, true, 4 )
 482 {
 483 }
 484
 485
 486 void Core_DotProductTest::run_func()
 487 {
 488     test_mat[OUTPUT][0].at<Scalar>(0,0) = Scalar(cvDotProduct( test_array[INPUT][0], test_array[INPUT][1] ));
 489 }
 490
 491
 492 void Core_DotProductTest::prepare_to_validation( int )
 493 {
 494     test_mat[REF_OUTPUT][0].at<Scalar>(0,0) = Scalar(cvtest::crossCorr( test_mat[INPUT][0], test_mat[INPUT][1] ));
 495 }
 496
 497
 498 ///////// crossproduct //////////
 499
 500 class Core_CrossProductTest : public Core_MatrixTest
 501 {
 502 public:
 503     Core_CrossProductTest();
 504 protected:
 505     void get_test_array_types_and_sizes( int test_case_idx,
 506                                         vector<vector<Size> >& sizes,
 507                                         vector<vector<int> >& types );
 508     void run_func();
 509     void prepare_to_validation( int test_case_idx );
 510 };
 511
 512
 513 Core_CrossProductTest::Core_CrossProductTest() : Core_MatrixTest( 2, 1, false, false, 1 )
 514 {
 515 }
 516
 517
 518 void Core_CrossProductTest::get_test_array_types_and_sizes( int,
 519                                                              vector<vector<Size> >& sizes,
 520                                                              vector<vector<int> >& types )
 521 {
 522     RNG& rng = ts->get_rng();
 523     int depth = cvtest::randInt(rng) % 2 + CV_32F;
 524     int cn = cvtest::randInt(rng) & 1 ? 3 : 1, type = CV_MAKETYPE(depth, cn);
 525     CvSize sz;
 526
 527     types[INPUT][0] = types[INPUT][1] = types[OUTPUT][0] = types[REF_OUTPUT][0] = type;
 528
 529     if( cn == 3 )
 530         sz = Size(1,1);
 531     else if( cvtest::randInt(rng) & 1 )
 532         sz = Size(3,1);
 533     else
 534         sz = Size(1,3);
 535
 536     sizes[INPUT][0] = sizes[INPUT][1] = sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = sz;
 537 }
 538
 539
 540 void Core_CrossProductTest::run_func()
 541 {
 542     cvCrossProduct( test_array[INPUT][0], test_array[INPUT][1], test_array[OUTPUT][0] );
 543 }
 544
 545
 546 void Core_CrossProductTest::prepare_to_validation( int )
 547 {
 548     CvScalar a = {{0,0,0,0}}, b = {{0,0,0,0}}, c = {{0,0,0,0}};
 549
 550     if( test_mat[INPUT][0].rows > 1 )
 551     {
 552         a.val[0] = cvGetReal2D( test_array[INPUT][0], 0, 0 );
 553         a.val[1] = cvGetReal2D( test_array[INPUT][0], 1, 0 );
 554         a.val[2] = cvGetReal2D( test_array[INPUT][0], 2, 0 );
 555
 556         b.val[0] = cvGetReal2D( test_array[INPUT][1], 0, 0 );
 557         b.val[1] = cvGetReal2D( test_array[INPUT][1], 1, 0 );
 558         b.val[2] = cvGetReal2D( test_array[INPUT][1], 2, 0 );
 559     }
 560     else if( test_mat[INPUT][0].cols > 1 )
 561     {
 562         a.val[0] = cvGetReal1D( test_array[INPUT][0], 0 );
 563         a.val[1] = cvGetReal1D( test_array[INPUT][0], 1 );
 564         a.val[2] = cvGetReal1D( test_array[INPUT][0], 2 );
 565
 566         b.val[0] = cvGetReal1D( test_array[INPUT][1], 0 );
 567         b.val[1] = cvGetReal1D( test_array[INPUT][1], 1 );
 568         b.val[2] = cvGetReal1D( test_array[INPUT][1], 2 );
 569     }
 570     else
 571     {
 572         a = cvGet1D( test_array[INPUT][0], 0 );
 573         b = cvGet1D( test_array[INPUT][1], 0 );
 574     }
 575
 576     c.val[2] = a.val[0]*b.val[1] - a.val[1]*b.val[0];
 577     c.val[1] = -a.val[0]*b.val[2] + a.val[2]*b.val[0];
 578     c.val[0] = a.val[1]*b.val[2] - a.val[2]*b.val[1];
 579
 580     if( test_mat[REF_OUTPUT][0].rows > 1 )
 581     {
 582         cvSetReal2D( test_array[REF_OUTPUT][0], 0, 0, c.val[0] );
 583         cvSetReal2D( test_array[REF_OUTPUT][0], 1, 0, c.val[1] );
 584         cvSetReal2D( test_array[REF_OUTPUT][0], 2, 0, c.val[2] );
 585     }
 586     else if( test_mat[REF_OUTPUT][0].cols > 1 )
 587     {
 588         cvSetReal1D( test_array[REF_OUTPUT][0], 0, c.val[0] );
 589         cvSetReal1D( test_array[REF_OUTPUT][0], 1, c.val[1] );
 590         cvSetReal1D( test_array[REF_OUTPUT][0], 2, c.val[2] );
 591     }
 592     else
 593     {
 594         cvSet1D( test_array[REF_OUTPUT][0], 0, c );
 595     }
 596 }
 597
 598
 599 ///////////////// gemm /////////////////////
 600
 601 class Core_GEMMTest : public Core_MatrixTest
 602 {
 603 public:
 604     typedef Core_MatrixTest Base;
 605     Core_GEMMTest();
 606 protected:
 607     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
 608     void get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high );
 609     int prepare_test_case( int test_case_idx );
 610     void run_func();
 611     void prepare_to_validation( int test_case_idx );
 612     int tabc_flag;
 613     double alpha, beta;
 614 };
 615
 616 Core_GEMMTest::Core_GEMMTest() : Core_MatrixTest( 5, 1, false, false, 2 )
 617 {
 618     test_case_count = 100;
 619     max_log_array_size = 10;
 620     alpha = beta = 0;
 621 }
 622
 623
 624 void Core_GEMMTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
 625 {
 626     RNG& rng = ts->get_rng();
 627     Size sizeA;
 628     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
 629     sizeA = sizes[INPUT][0];
 630     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
 631     sizes[INPUT][0] = sizeA;
 632     sizes[INPUT][2] = sizes[INPUT][3] = Size(1,1);
 633     types[INPUT][2] = types[INPUT][3] &= ~CV_MAT_CN_MASK;
 634
 635     tabc_flag = cvtest::randInt(rng) & 7;
 636
 637     switch( tabc_flag & (CV_GEMM_A_T|CV_GEMM_B_T) )
 638     {
 639         case 0:
 640             sizes[INPUT][1].height = sizes[INPUT][0].width;
 641             sizes[OUTPUT][0].height = sizes[INPUT][0].height;
 642             sizes[OUTPUT][0].width = sizes[INPUT][1].width;
 643             break;
 644         case CV_GEMM_B_T:
 645             sizes[INPUT][1].width = sizes[INPUT][0].width;
 646             sizes[OUTPUT][0].height = sizes[INPUT][0].height;
 647             sizes[OUTPUT][0].width = sizes[INPUT][1].height;
 648             break;
 649         case CV_GEMM_A_T:
 650             sizes[INPUT][1].height = sizes[INPUT][0].height;
 651             sizes[OUTPUT][0].height = sizes[INPUT][0].width;
 652             sizes[OUTPUT][0].width = sizes[INPUT][1].width;
 653             break;
 654         case CV_GEMM_A_T | CV_GEMM_B_T:
 655             sizes[INPUT][1].width = sizes[INPUT][0].height;
 656             sizes[OUTPUT][0].height = sizes[INPUT][0].width;
 657             sizes[OUTPUT][0].width = sizes[INPUT][1].height;
 658             break;
 659     }
 660
 661     sizes[REF_OUTPUT][0] = sizes[OUTPUT][0];
 662
 663     if( cvtest::randInt(rng) & 1 )
 664         sizes[INPUT][4] = Size(0,0);
 665     else if( !(tabc_flag & CV_GEMM_C_T) )
 666         sizes[INPUT][4] = sizes[OUTPUT][0];
 667     else
 668     {
 669         sizes[INPUT][4].width = sizes[OUTPUT][0].height;
 670         sizes[INPUT][4].height = sizes[OUTPUT][0].width;
 671     }
 672 }
 673
 674
 675 int Core_GEMMTest::prepare_test_case( int test_case_idx )
 676 {
 677     int code = Base::prepare_test_case( test_case_idx );
 678     if( code > 0 )
 679     {
 680         alpha = cvGetReal2D( test_array[INPUT][2], 0, 0 );
 681         beta = cvGetReal2D( test_array[INPUT][3], 0, 0 );
 682     }
 683     return code;
 684 }
 685
 686
 687 void Core_GEMMTest::get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high )
 688 {
 689     low = Scalar::all(-10.);
 690     high = Scalar::all(10.);
 691 }
 692
 693
 694 void Core_GEMMTest::run_func()
 695 {
 696     cvGEMM( test_array[INPUT][0], test_array[INPUT][1], alpha,
 697            test_array[INPUT][4], beta, test_array[OUTPUT][0], tabc_flag );
 698 }
 699
 700
 701 void Core_GEMMTest::prepare_to_validation( int )
 702 {
 703     cvtest::gemm( test_mat[INPUT][0], test_mat[INPUT][1], alpha,
 704              test_array[INPUT][4] ? test_mat[INPUT][4] : Mat(),
 705              beta, test_mat[REF_OUTPUT][0], tabc_flag );
 706 }
 707
 708
 709 ///////////////// multransposed /////////////////////
 710
 711 class Core_MulTransposedTest : public Core_MatrixTest
 712 {
 713 public:
 714     Core_MulTransposedTest();
 715 protected:
 716     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
 717     void get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high );
 718     void run_func();
 719     void prepare_to_validation( int test_case_idx );
 720     int order;
 721 };
 722
 723
 724 Core_MulTransposedTest::Core_MulTransposedTest() : Core_MatrixTest( 2, 1, false, false, 1 )
 725 {
 726     test_case_count = 100;
 727     order = 0;
 728     test_array[TEMP].push_back(NULL);
 729 }
 730
 731
 732 void Core_MulTransposedTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
 733 {
 734     RNG& rng = ts->get_rng();
 735     int bits = cvtest::randInt(rng);
 736     int src_type = cvtest::randInt(rng) % 5;
 737     int dst_type = cvtest::randInt(rng) % 2;
 738
 739     src_type = src_type == 0 ? CV_8U : src_type == 1 ? CV_16U : src_type == 2 ? CV_16S :
 740     src_type == 3 ? CV_32F : CV_64F;
 741     dst_type = dst_type == 0 ? CV_32F : CV_64F;
 742     dst_type = MAX( dst_type, src_type );
 743
 744     Core_MatrixTest::get_test_array_types_and_sizes( test_case_idx, sizes, types );
 745
 746     if( bits & 1 )
 747         sizes[INPUT][1] = Size(0,0);
 748     else
 749     {
 750         sizes[INPUT][1] = sizes[INPUT][0];
 751         if( bits & 2 )
 752             sizes[INPUT][1].height = 1;
 753         if( bits & 4 )
 754             sizes[INPUT][1].width = 1;
 755     }
 756
 757     sizes[TEMP][0] = sizes[INPUT][0];
 758     types[INPUT][0] = src_type;
 759     types[OUTPUT][0] = types[REF_OUTPUT][0] = types[INPUT][1] = types[TEMP][0] = dst_type;
 760
 761     order = (bits & 8) != 0;
 762     sizes[OUTPUT][0].width = sizes[OUTPUT][0].height = order == 0 ?
 763     sizes[INPUT][0].height : sizes[INPUT][0].width;
 764     sizes[REF_OUTPUT][0] = sizes[OUTPUT][0];
 765 }
 766
 767
 768 void Core_MulTransposedTest::get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high )
 769 {
 770     low = cvScalarAll(-10.);
 771     high = cvScalarAll(10.);
 772 }
 773
 774
 775 void Core_MulTransposedTest::run_func()
 776 {
 777     cvMulTransposed( test_array[INPUT][0], test_array[OUTPUT][0],
 778                     order, test_array[INPUT][1] );
 779 }
 780
 781
 782 void Core_MulTransposedTest::prepare_to_validation( int )
 783 {
 784     const Mat& src = test_mat[INPUT][0];
 785     Mat delta = test_mat[INPUT][1];
 786     Mat& temp = test_mat[TEMP][0];
 787     if( !delta.empty() )
 788     {
 789         if( delta.rows < src.rows || delta.cols < src.cols )
 790         {
 791             cv::repeat( delta, src.rows/delta.rows, src.cols/delta.cols, temp);
 792             delta = temp;
 793         }
 794         cvtest::add( src, 1, delta, -1, Scalar::all(0), temp, temp.type());
 795     }
 796     else
 797         src.convertTo(temp, temp.type());
 798
 799     cvtest::gemm( temp, temp, 1., Mat(), 0, test_mat[REF_OUTPUT][0], order == 0 ? GEMM_2_T : GEMM_1_T );
 800 }
 801
 802
 803 ///////////////// Transform /////////////////////
 804
 805 class Core_TransformTest : public Core_MatrixTest
 806 {
 807 public:
 808     typedef Core_MatrixTest Base;
 809     Core_TransformTest();
 810 protected:
 811     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
 812     double get_success_error_level( int test_case_idx, int i, int j );
 813     int prepare_test_case( int test_case_idx );
 814     void run_func();
 815     void prepare_to_validation( int test_case_idx );
 816
 817     double scale;
 818     bool diagMtx;
 819 };
 820
 821
 822 Core_TransformTest::Core_TransformTest() : Core_MatrixTest( 3, 1, true, false, 4 )
 823 {
 824 }
 825
 826
 827 void Core_TransformTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
 828 {
 829     RNG& rng = ts->get_rng();
 830     int bits = cvtest::randInt(rng);
 831     int depth, dst_cn, mat_cols, mattype;
 832     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
 833
 834     mat_cols = CV_MAT_CN(types[INPUT][0]);
 835     depth = CV_MAT_DEPTH(types[INPUT][0]);
 836     dst_cn = cvtest::randInt(rng) % 4 + 1;
 837     types[OUTPUT][0] = types[REF_OUTPUT][0] = CV_MAKETYPE(depth, dst_cn);
 838
 839     mattype = depth < CV_32S ? CV_32F : depth == CV_64F ? CV_64F : bits & 1 ? CV_32F : CV_64F;
 840     types[INPUT][1] = mattype;
 841     types[INPUT][2] = CV_MAKETYPE(mattype, dst_cn);
 842
 843     scale = 1./((cvtest::randInt(rng)%4)*50+1);
 844
 845     if( bits & 2 )
 846     {
 847         sizes[INPUT][2] = Size(0,0);
 848         mat_cols += (bits & 4) != 0;
 849     }
 850     else if( bits & 4 )
 851         sizes[INPUT][2] = Size(1,1);
 852     else
 853     {
 854         if( bits & 8 )
 855             sizes[INPUT][2] = Size(dst_cn,1);
 856         else
 857             sizes[INPUT][2] = Size(1,dst_cn);
 858         types[INPUT][2] &= ~CV_MAT_CN_MASK;
 859     }
 860     diagMtx = (bits & 16) != 0;
 861
 862     sizes[INPUT][1] = Size(mat_cols,dst_cn);
 863 }
 864
 865
 866 int Core_TransformTest::prepare_test_case( int test_case_idx )
 867 {
 868     int code = Base::prepare_test_case( test_case_idx );
 869     if( code > 0 )
 870     {
 871         Mat& m = test_mat[INPUT][1];
 872         cvtest::add(m, scale, m, 0, Scalar::all(0), m, m.type() );
 873         if(diagMtx)
 874         {
 875             Mat mask = Mat::eye(m.rows, m.cols, CV_8U)*255;
 876             mask = ~mask;
 877             m.setTo(Scalar::all(0), mask);
 878         }
 879     }
 880     return code;
 881 }
 882
 883
 884 double Core_TransformTest::get_success_error_level( int test_case_idx, int i, int j )
 885 {
 886     int depth = test_mat[INPUT][0].depth();
 887     return depth <= CV_8S ? 1 : depth <= CV_32S ? 9 : Base::get_success_error_level( test_case_idx, i, j );
 888 }
 889
 890 void Core_TransformTest::run_func()
 891 {
 892     CvMat _m = test_mat[INPUT][1], _shift = test_mat[INPUT][2];
 893     cvTransform( test_array[INPUT][0], test_array[OUTPUT][0], &_m, _shift.data.ptr ? &_shift : 0);
 894 }
 895
 896
 897 void Core_TransformTest::prepare_to_validation( int )
 898 {
 899     Mat transmat = test_mat[INPUT][1];
 900     Mat shift = test_mat[INPUT][2];
 901
 902     cvtest::transform( test_mat[INPUT][0], test_mat[REF_OUTPUT][0], transmat, shift );
 903 }
 904
 905
 906 ///////////////// PerspectiveTransform /////////////////////
 907
 908 class Core_PerspectiveTransformTest : public Core_MatrixTest
 909 {
 910 public:
 911     Core_PerspectiveTransformTest();
 912 protected:
 913     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
 914     double get_success_error_level( int test_case_idx, int i, int j );
 915     void run_func();
 916     void prepare_to_validation( int test_case_idx );
 917 };
 918
 919
 920 Core_PerspectiveTransformTest::Core_PerspectiveTransformTest() : Core_MatrixTest( 2, 1, false, false, 2 )
 921 {
 922 }
 923
 924
 925 void Core_PerspectiveTransformTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
 926 {
 927     RNG& rng = ts->get_rng();
 928     int bits = cvtest::randInt(rng);
 929     int depth, cn, mattype;
 930     Core_MatrixTest::get_test_array_types_and_sizes( test_case_idx, sizes, types );
 931
 932     cn = CV_MAT_CN(types[INPUT][0]) + 1;
 933     depth = CV_MAT_DEPTH(types[INPUT][0]);
 934     types[INPUT][0] = types[OUTPUT][0] = types[REF_OUTPUT][0] = CV_MAKETYPE(depth, cn);
 935
 936     mattype = depth == CV_64F ? CV_64F : bits & 1 ? CV_32F : CV_64F;
 937     types[INPUT][1] = mattype;
 938     sizes[INPUT][1] = Size(cn + 1, cn + 1);
 939 }
 940
 941
 942 double Core_PerspectiveTransformTest::get_success_error_level( int test_case_idx, int i, int j )
 943 {
 944     int depth = test_mat[INPUT][0].depth();
 945     return depth == CV_32F ? 1e-4 : depth == CV_64F ? 1e-8 :
 946     Core_MatrixTest::get_success_error_level(test_case_idx, i, j);
 947 }
 948
 949
 950 void Core_PerspectiveTransformTest::run_func()
 951 {
 952     CvMat _m = test_mat[INPUT][1];
 953     cvPerspectiveTransform( test_array[INPUT][0], test_array[OUTPUT][0], &_m );
 954 }
 955
 956
 957 static void cvTsPerspectiveTransform( const CvArr* _src, CvArr* _dst, const CvMat* transmat )
 958 {
 959     int i, j, cols;
 960     int cn, depth, mat_depth;
 961     CvMat astub, bstub, *a, *b;
 962     double mat[16];
 963
 964     a = cvGetMat( _src, &astub, 0, 0 );
 965     b = cvGetMat( _dst, &bstub, 0, 0 );
 966
 967     cn = CV_MAT_CN(a->type);
 968     depth = CV_MAT_DEPTH(a->type);
 969     mat_depth = CV_MAT_DEPTH(transmat->type);
 970     cols = transmat->cols;
 971
 972     // prepare cn x (cn + 1) transform matrix
 973     if( mat_depth == CV_32F )
 974     {
 975         for( i = 0; i < transmat->rows; i++ )
 976             for( j = 0; j < cols; j++ )
 977                 mat[i*cols + j] = ((float*)(transmat->data.ptr + transmat->step*i))[j];
 978     }
 979     else
 980     {
 981         assert( mat_depth == CV_64F );
 982         for( i = 0; i < transmat->rows; i++ )
 983             for( j = 0; j < cols; j++ )
 984                 mat[i*cols + j] = ((double*)(transmat->data.ptr + transmat->step*i))[j];
 985     }
 986
 987     // transform data
 988     cols = a->cols * cn;
 989     vector<double> buf(cols);
 990
 991     for( i = 0; i < a->rows; i++ )
 992     {
 993         uchar* src = a->data.ptr + i*a->step;
 994         uchar* dst = b->data.ptr + i*b->step;
 995
 996         switch( depth )
 997         {
 998             case CV_32F:
 999                 for( j = 0; j < cols; j++ )
1000                     buf[j] = ((float*)src)[j];
1001                 break;
1002             case CV_64F:
1003                 for( j = 0; j < cols; j++ )
1004                     buf[j] = ((double*)src)[j];
1005                 break;
1006             default:
1007                 assert(0);
1008         }
1009
1010         switch( cn )
1011         {
1012             case 2:
1013                 for( j = 0; j < cols; j += 2 )
1014                 {
1015                     double t0 = buf[j]*mat[0] + buf[j+1]*mat[1] + mat[2];
1016                     double t1 = buf[j]*mat[3] + buf[j+1]*mat[4] + mat[5];
1017                     double w = buf[j]*mat[6] + buf[j+1]*mat[7] + mat[8];
1018                     w = w ? 1./w : 0;
1019                     buf[j] = t0*w;
1020                     buf[j+1] = t1*w;
1021                 }
1022                 break;
1023             case 3:
1024                 for( j = 0; j < cols; j += 3 )
1025                 {
1026                     double t0 = buf[j]*mat[0] + buf[j+1]*mat[1] + buf[j+2]*mat[2] + mat[3];
1027                     double t1 = buf[j]*mat[4] + buf[j+1]*mat[5] + buf[j+2]*mat[6] + mat[7];
1028                     double t2 = buf[j]*mat[8] + buf[j+1]*mat[9] + buf[j+2]*mat[10] + mat[11];
1029                     double w = buf[j]*mat[12] + buf[j+1]*mat[13] + buf[j+2]*mat[14] + mat[15];
1030                     w = w ? 1./w : 0;
1031                     buf[j] = t0*w;
1032                     buf[j+1] = t1*w;
1033                     buf[j+2] = t2*w;
1034                 }
1035                 break;
1036             default:
1037                 assert(0);
1038         }
1039
1040         switch( depth )
1041         {
1042             case CV_32F:
1043                 for( j = 0; j < cols; j++ )
1044                     ((float*)dst)[j] = (float)buf[j];
1045                 break;
1046             case CV_64F:
1047                 for( j = 0; j < cols; j++ )
1048                     ((double*)dst)[j] = buf[j];
1049                 break;
1050             default:
1051                 assert(0);
1052         }
1053     }
1054 }
1055
1056
1057 void Core_PerspectiveTransformTest::prepare_to_validation( int )
1058 {
1059     CvMat transmat = test_mat[INPUT][1];
1060     cvTsPerspectiveTransform( test_array[INPUT][0], test_array[REF_OUTPUT][0], &transmat );
1061 }
1062
1063 ///////////////// Mahalanobis /////////////////////
1064
1065 class Core_MahalanobisTest : public Core_MatrixTest
1066 {
1067 public:
1068     typedef Core_MatrixTest Base;
1069     Core_MahalanobisTest();
1070 protected:
1071     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1072     int prepare_test_case( int test_case_idx );
1073     void run_func();
1074     void prepare_to_validation( int test_case_idx );
1075 };
1076
1077
1078 Core_MahalanobisTest::Core_MahalanobisTest() : Core_MatrixTest( 3, 1, false, true, 1 )
1079 {
1080     test_case_count = 100;
1081     test_array[TEMP].push_back(NULL);
1082     test_array[TEMP].push_back(NULL);
1083     test_array[TEMP].push_back(NULL);
1084 }
1085
1086
1087 void Core_MahalanobisTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
1088 {
1089     RNG& rng = ts->get_rng();
1090     Core_MatrixTest::get_test_array_types_and_sizes( test_case_idx, sizes, types );
1091
1092     if( cvtest::randInt(rng) & 1 )
1093         sizes[INPUT][0].width = sizes[INPUT][1].width = 1;
1094     else
1095         sizes[INPUT][0].height = sizes[INPUT][1].height = 1;
1096
1097     sizes[TEMP][0] = sizes[TEMP][1] = sizes[INPUT][0];
1098     sizes[INPUT][2].width = sizes[INPUT][2].height = sizes[INPUT][0].width + sizes[INPUT][0].height - 1;
1099     sizes[TEMP][2] = sizes[INPUT][2];
1100     types[TEMP][0] = types[TEMP][1] = types[TEMP][2] = types[INPUT][0];
1101 }
1102
1103 int Core_MahalanobisTest::prepare_test_case( int test_case_idx )
1104 {
1105     int code = Base::prepare_test_case( test_case_idx );
1106     if( code > 0 )
1107     {
1108         // make sure that the inverted "covariation" matrix is symmetrix and positively defined.
1109         cvtest::gemm( test_mat[INPUT][2], test_mat[INPUT][2], 1., Mat(), 0., test_mat[TEMP][2], GEMM_2_T );
1110         cvtest::copy( test_mat[TEMP][2], test_mat[INPUT][2] );
1111     }
1112
1113     return code;
1114 }
1115
1116
1117 void Core_MahalanobisTest::run_func()
1118 {
1119     test_mat[OUTPUT][0].at<Scalar>(0,0) =
1120     cvRealScalar(cvMahalanobis(test_array[INPUT][0], test_array[INPUT][1], test_array[INPUT][2]));
1121 }
1122
1123 void Core_MahalanobisTest::prepare_to_validation( int )
1124 {
1125     cvtest::add( test_mat[INPUT][0], 1., test_mat[INPUT][1], -1.,
1126                 Scalar::all(0), test_mat[TEMP][0], test_mat[TEMP][0].type() );
1127     if( test_mat[INPUT][0].rows == 1 )
1128         cvtest::gemm( test_mat[TEMP][0], test_mat[INPUT][2], 1.,
1129                  Mat(), 0., test_mat[TEMP][1], 0 );
1130     else
1131         cvtest::gemm( test_mat[INPUT][2], test_mat[TEMP][0], 1.,
1132                  Mat(), 0., test_mat[TEMP][1], 0 );
1133
1134     test_mat[REF_OUTPUT][0].at<Scalar>(0,0) = cvRealScalar(sqrt(cvtest::crossCorr(test_mat[TEMP][0], test_mat[TEMP][1])));
1135 }
1136
1137
1138 ///////////////// covarmatrix /////////////////////
1139
1140 class Core_CovarMatrixTest : public Core_MatrixTest
1141 {
1142 public:
1143     Core_CovarMatrixTest();
1144 protected:
1145     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1146     int prepare_test_case( int test_case_idx );
1147     void run_func();
1148     void prepare_to_validation( int test_case_idx );
1149     vector<void*> temp_hdrs;
1150     vector<uchar> hdr_data;
1151     int flags, t_flag, len, count;
1152     bool are_images;
1153 };
1154
1155
1156 Core_CovarMatrixTest::Core_CovarMatrixTest() : Core_MatrixTest( 1, 1, true, false, 1 ),
1157 flags(0), t_flag(0), are_images(false)
1158 {
1159     test_case_count = 100;
1160     test_array[INPUT_OUTPUT].push_back(NULL);
1161     test_array[REF_INPUT_OUTPUT].push_back(NULL);
1162     test_array[TEMP].push_back(NULL);
1163     test_array[TEMP].push_back(NULL);
1164 }
1165
1166
1167 void Core_CovarMatrixTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
1168 {
1169     RNG& rng = ts->get_rng();
1170     int bits = cvtest::randInt(rng);
1171     int i, single_matrix;
1172     Core_MatrixTest::get_test_array_types_and_sizes( test_case_idx, sizes, types );
1173
1174     flags = bits & (CV_COVAR_NORMAL | CV_COVAR_USE_AVG | CV_COVAR_SCALE | CV_COVAR_ROWS );
1175     single_matrix = flags & CV_COVAR_ROWS;
1176     t_flag = (bits & 256) != 0;
1177
1178     const int min_count = 2;
1179
1180     if( !t_flag )
1181     {
1182         len = sizes[INPUT][0].width;
1183         count = sizes[INPUT][0].height;
1184         count = MAX(count, min_count);
1185         sizes[INPUT][0] = Size(len, count);
1186     }
1187     else
1188     {
1189         len = sizes[INPUT][0].height;
1190         count = sizes[INPUT][0].width;
1191         count = MAX(count, min_count);
1192         sizes[INPUT][0] = Size(count, len);
1193     }
1194
1195     if( single_matrix && t_flag )
1196         flags = (flags & ~CV_COVAR_ROWS) | CV_COVAR_COLS;
1197
1198     if( CV_MAT_DEPTH(types[INPUT][0]) == CV_32S )
1199         types[INPUT][0] = (types[INPUT][0] & ~CV_MAT_DEPTH_MASK) | CV_32F;
1200
1201     sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = flags & CV_COVAR_NORMAL ? Size(len,len) : Size(count,count);
1202     sizes[INPUT_OUTPUT][0] = sizes[REF_INPUT_OUTPUT][0] = !t_flag ? Size(len,1) : Size(1,len);
1203     sizes[TEMP][0] = sizes[INPUT][0];
1204
1205     types[INPUT_OUTPUT][0] = types[REF_INPUT_OUTPUT][0] =
1206     types[OUTPUT][0] = types[REF_OUTPUT][0] = types[TEMP][0] =
1207     CV_MAT_DEPTH(types[INPUT][0]) == CV_64F || (bits & 512) ? CV_64F : CV_32F;
1208
1209     are_images = (bits & 1024) != 0;
1210     for( i = 0; i < (single_matrix ? 1 : count); i++ )
1211         temp_hdrs.push_back(NULL);
1212 }
1213
1214
1215 int Core_CovarMatrixTest::prepare_test_case( int test_case_idx )
1216 {
1217     int code = Core_MatrixTest::prepare_test_case( test_case_idx );
1218     if( code > 0 )
1219     {
1220         int i;
1221         int single_matrix = flags & (CV_COVAR_ROWS|CV_COVAR_COLS);
1222         int hdr_size = are_images ? sizeof(IplImage) : sizeof(CvMat);
1223
1224         hdr_data.resize(count*hdr_size);
1225         uchar* _hdr_data = &hdr_data[0];
1226         if( single_matrix )
1227         {
1228             if( !are_images )
1229                 *((CvMat*)_hdr_data) = test_mat[INPUT][0];
1230             else
1231                 *((IplImage*)_hdr_data) = test_mat[INPUT][0];
1232             temp_hdrs[0] = _hdr_data;
1233         }
1234         else
1235             for( i = 0; i < count; i++ )
1236             {
1237                 Mat part;
1238                 void* ptr = _hdr_data + i*hdr_size;
1239
1240                 if( !t_flag )
1241                     part = test_mat[INPUT][0].row(i);
1242                 else
1243                     part = test_mat[INPUT][0].col(i);
1244
1245                 if( !are_images )
1246                     *((CvMat*)ptr) = part;
1247                 else
1248                     *((IplImage*)ptr) = part;
1249
1250                 temp_hdrs[i] = ptr;
1251             }
1252     }
1253
1254     return code;
1255 }
1256
1257
1258 void Core_CovarMatrixTest::run_func()
1259 {
1260     cvCalcCovarMatrix( (const void**)&temp_hdrs[0], count,
1261                       test_array[OUTPUT][0], test_array[INPUT_OUTPUT][0], flags );
1262 }
1263
1264
1265 void Core_CovarMatrixTest::prepare_to_validation( int )
1266 {
1267     Mat& avg = test_mat[REF_INPUT_OUTPUT][0];
1268     double scale = 1.;
1269
1270     if( !(flags & CV_COVAR_USE_AVG) )
1271     {
1272         Mat hdrs0 = cvarrToMat(temp_hdrs[0]);
1273
1274         int i;
1275         avg = Scalar::all(0);
1276
1277         for( i = 0; i < count; i++ )
1278         {
1279             Mat vec;
1280             if( flags & CV_COVAR_ROWS )
1281                 vec = hdrs0.row(i);
1282             else if( flags & CV_COVAR_COLS )
1283                 vec = hdrs0.col(i);
1284             else
1285                 vec = cvarrToMat(temp_hdrs[i]);
1286
1287             cvtest::add(avg, 1, vec, 1, Scalar::all(0), avg, avg.type());
1288         }
1289
1290         cvtest::add(avg, 1./count, avg, 0., Scalar::all(0), avg, avg.type());
1291     }
1292
1293     if( flags & CV_COVAR_SCALE )
1294     {
1295         scale = 1./count;
1296     }
1297
1298     Mat& temp0 = test_mat[TEMP][0];
1299     cv::repeat( avg, temp0.rows/avg.rows, temp0.cols/avg.cols, temp0 );
1300     cvtest::add( test_mat[INPUT][0], 1, temp0, -1, Scalar::all(0), temp0, temp0.type());
1301
1302     cvtest::gemm( temp0, temp0, scale, Mat(), 0., test_mat[REF_OUTPUT][0],
1303              t_flag ^ ((flags & CV_COVAR_NORMAL) != 0) ? CV_GEMM_A_T : CV_GEMM_B_T );
1304     temp_hdrs.clear();
1305 }
1306
1307
1308 static void cvTsFloodWithZeros( Mat& mat, RNG& rng )
1309 {
1310     int k, total = mat.rows*mat.cols, type = mat.type();
1311     int zero_total = cvtest::randInt(rng) % total;
1312     CV_Assert( type == CV_32FC1 || type == CV_64FC1 );
1313
1314     for( k = 0; k < zero_total; k++ )
1315     {
1316         int i = cvtest::randInt(rng) % mat.rows;
1317         int j = cvtest::randInt(rng) % mat.cols;
1318
1319         if( type == CV_32FC1 )
1320             mat.at<float>(i,j) = 0.f;
1321         else
1322             mat.at<double>(i,j) = 0.;
1323     }
1324 }
1325
1326
1327 ///////////////// determinant /////////////////////
1328
1329 class Core_DetTest : public Core_MatrixTest
1330 {
1331 public:
1332     typedef Core_MatrixTest Base;
1333     Core_DetTest();
1334 protected:
1335     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1336     double get_success_error_level( int test_case_idx, int i, int j );
1337     void get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high );
1338     int prepare_test_case( int test_case_idx );
1339     void run_func();
1340     void prepare_to_validation( int test_case_idx );
1341 };
1342
1343
1344 Core_DetTest::Core_DetTest() : Core_MatrixTest( 1, 1, false, true, 1 )
1345 {
1346     test_case_count = 100;
1347     max_log_array_size = 7;
1348     test_array[TEMP].push_back(NULL);
1349 }
1350
1351
1352 void Core_DetTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
1353 {
1354     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
1355
1356     sizes[INPUT][0].width = sizes[INPUT][0].height = sizes[INPUT][0].height;
1357     sizes[TEMP][0] = sizes[INPUT][0];
1358     types[TEMP][0] = CV_64FC1;
1359 }
1360
1361
1362 void Core_DetTest::get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high )
1363 {
1364     low = cvScalarAll(-2.);
1365     high = cvScalarAll(2.);
1366 }
1367
1368
1369 double Core_DetTest::get_success_error_level( int /*test_case_idx*/, int /*i*/, int /*j*/ )
1370 {
1371     return CV_MAT_DEPTH(cvGetElemType(test_array[INPUT][0])) == CV_32F ? 1e-2 : 1e-5;
1372 }
1373
1374
1375 int Core_DetTest::prepare_test_case( int test_case_idx )
1376 {
1377     int code = Core_MatrixTest::prepare_test_case( test_case_idx );
1378     if( code > 0 )
1379         cvTsFloodWithZeros( test_mat[INPUT][0], ts->get_rng() );
1380
1381     return code;
1382 }
1383
1384
1385 void Core_DetTest::run_func()
1386 {
1387     test_mat[OUTPUT][0].at<Scalar>(0,0) = cvRealScalar(cvDet(test_array[INPUT][0]));
1388 }
1389
1390
1391 // LU method that chooses the optimal in a column pivot element
1392 static double cvTsLU( CvMat* a, CvMat* b=NULL, CvMat* x=NULL, int* rank=0 )
1393 {
1394     int i, j, k, N = a->rows, N1 = a->cols, Nm = MIN(N, N1), step = a->step/sizeof(double);
1395     int M = b ? b->cols : 0, b_step = b ? b->step/sizeof(double) : 0;
1396     int x_step = x ? x->step/sizeof(double) : 0;
1397     double *a0 = a->data.db, *b0 = b ? b->data.db : 0;
1398     double *x0 = x ? x->data.db : 0;
1399     double t, det = 1.;
1400     assert( CV_MAT_TYPE(a->type) == CV_64FC1 &&
1401            (!b || CV_ARE_TYPES_EQ(a,b)) && (!x || CV_ARE_TYPES_EQ(a,x)));
1402
1403     for( i = 0; i < Nm; i++ )
1404     {
1405         double max_val = fabs(a0[i*step + i]);
1406         double *a1, *a2, *b1 = 0, *b2 = 0;
1407         k = i;
1408
1409         for( j = i+1; j < N; j++ )
1410         {
1411             t = fabs(a0[j*step + i]);
1412             if( max_val < t )
1413             {
1414                 max_val = t;
1415                 k = j;
1416             }
1417         }
1418
1419         if( k != i )
1420         {
1421             for( j = i; j < N1; j++ )
1422                 CV_SWAP( a0[i*step + j], a0[k*step + j], t );
1423
1424             for( j = 0; j < M; j++ )
1425                 CV_SWAP( b0[i*b_step + j], b0[k*b_step + j], t );
1426             det = -det;
1427         }
1428
1429         if( max_val == 0 )
1430         {
1431             if( rank )
1432                 *rank = i;
1433             return 0.;
1434         }
1435
1436         a1 = a0 + i*step;
1437         a2 = a1 + step;
1438         b1 = b0 + i*b_step;
1439         b2 = b1 + b_step;
1440
1441         for( j = i+1; j < N; j++, a2 += step, b2 += b_step )
1442         {
1443             t = a2[i]/a1[i];
1444             for( k = i+1; k < N1; k++ )
1445                 a2[k] -= t*a1[k];
1446
1447             for( k = 0; k < M; k++ )
1448                 b2[k] -= t*b1[k];
1449         }
1450
1451         det *= a1[i];
1452     }
1453
1454     if( x )
1455     {
1456         assert( b );
1457
1458         for( i = N-1; i >= 0; i-- )
1459         {
1460             double* a1 = a0 + i*step;
1461             double* b1 = b0 + i*b_step;
1462             for( j = 0; j < M; j++ )
1463             {
1464                 t = b1[j];
1465                 for( k = i+1; k < N1; k++ )
1466                     t -= a1[k]*x0[k*x_step + j];
1467                 x0[i*x_step + j] = t/a1[i];
1468             }
1469         }
1470     }
1471
1472     if( rank )
1473         *rank = i;
1474     return det;
1475 }
1476
1477
1478 void Core_DetTest::prepare_to_validation( int )
1479 {
1480     test_mat[INPUT][0].convertTo(test_mat[TEMP][0], test_mat[TEMP][0].type());
1481     CvMat temp0 = test_mat[TEMP][0];
1482     test_mat[REF_OUTPUT][0].at<Scalar>(0,0) = cvRealScalar(cvTsLU(&temp0, 0, 0));
1483 }
1484
1485
1486 ///////////////// invert /////////////////////
1487
1488 class Core_InvertTest : public Core_MatrixTest
1489 {
1490 public:
1491     typedef Core_MatrixTest Base;
1492     Core_InvertTest();
1493 protected:
1494     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1495     void get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high );
1496     double get_success_error_level( int test_case_idx, int i, int j );
1497     int prepare_test_case( int test_case_idx );
1498     void run_func();
1499     void prepare_to_validation( int test_case_idx );
1500     int method, rank;
1501     double result;
1502 };
1503
1504
1505 Core_InvertTest::Core_InvertTest()
1506 : Core_MatrixTest( 1, 1, false, false, 1 ), method(0), rank(0), result(0.)
1507 {
1508     test_case_count = 100;
1509     max_log_array_size = 7;
1510     test_array[TEMP].push_back(NULL);
1511     test_array[TEMP].push_back(NULL);
1512 }
1513
1514
1515 void Core_InvertTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
1516 {
1517     RNG& rng = ts->get_rng();
1518     int bits = cvtest::randInt(rng);
1519     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
1520     int min_size = MIN( sizes[INPUT][0].width, sizes[INPUT][0].height );
1521
1522     if( (bits & 3) == 0 )
1523     {
1524         method = CV_SVD;
1525         if( bits & 4 )
1526         {
1527             sizes[INPUT][0] = Size(min_size, min_size);
1528             if( bits & 16 )
1529                 method = CV_CHOLESKY;
1530         }
1531     }
1532     else
1533     {
1534         method = CV_LU;
1535         sizes[INPUT][0] = Size(min_size, min_size);
1536     }
1537
1538     sizes[TEMP][0].width = sizes[INPUT][0].height;
1539     sizes[TEMP][0].height = sizes[INPUT][0].width;
1540     sizes[TEMP][1] = sizes[INPUT][0];
1541     types[TEMP][0] = types[INPUT][0];
1542     types[TEMP][1] = CV_64FC1;
1543     sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = Size(min_size, min_size);
1544 }
1545
1546
1547 double Core_InvertTest::get_success_error_level( int /*test_case_idx*/, int, int )
1548 {
1549     return CV_MAT_DEPTH(cvGetElemType(test_array[OUTPUT][0])) == CV_32F ? 1e-2 : 1e-6;
1550 }
1551
1552 int Core_InvertTest::prepare_test_case( int test_case_idx )
1553 {
1554     int code = Core_MatrixTest::prepare_test_case( test_case_idx );
1555     if( code > 0 )
1556     {
1557         cvTsFloodWithZeros( test_mat[INPUT][0], ts->get_rng() );
1558
1559         if( method == CV_CHOLESKY )
1560         {
1561             cvtest::gemm( test_mat[INPUT][0], test_mat[INPUT][0], 1.,
1562                      Mat(), 0., test_mat[TEMP][0], CV_GEMM_B_T );
1563             cvtest::copy( test_mat[TEMP][0], test_mat[INPUT][0] );
1564         }
1565     }
1566
1567     return code;
1568 }
1569
1570
1571
1572 void Core_InvertTest::get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high )
1573 {
1574     low = cvScalarAll(-1.);
1575     high = cvScalarAll(1.);
1576 }
1577
1578
1579 void Core_InvertTest::run_func()
1580 {
1581     result = cvInvert(test_array[INPUT][0], test_array[TEMP][0], method);
1582 }
1583
1584
1585 static double cvTsSVDet( CvMat* mat, double* ratio )
1586 {
1587     int type = CV_MAT_TYPE(mat->type);
1588     int i, nm = MIN( mat->rows, mat->cols );
1589     CvMat* w = cvCreateMat( nm, 1, type );
1590     double det = 1.;
1591
1592     cvSVD( mat, w, 0, 0, 0 );
1593
1594     if( type == CV_32FC1 )
1595     {
1596         for( i = 0; i < nm; i++ )
1597             det *= w->data.fl[i];
1598         *ratio = w->data.fl[nm-1] < FLT_EPSILON ? 0 : w->data.fl[nm-1]/w->data.fl[0];
1599     }
1600     else
1601     {
1602         for( i = 0; i < nm; i++ )
1603             det *= w->data.db[i];
1604         *ratio = w->data.db[nm-1] < FLT_EPSILON ? 0 : w->data.db[nm-1]/w->data.db[0];
1605     }
1606
1607     cvReleaseMat( &w );
1608     return det;
1609 }
1610
1611 void Core_InvertTest::prepare_to_validation( int )
1612 {
1613     Mat& input = test_mat[INPUT][0];
1614     Mat& temp0 = test_mat[TEMP][0];
1615     Mat& temp1 = test_mat[TEMP][1];
1616     Mat& dst0 = test_mat[REF_OUTPUT][0];
1617     Mat& dst = test_mat[OUTPUT][0];
1618     CvMat _input = input;
1619     double ratio = 0, det = cvTsSVDet( &_input, &ratio );
1620     double threshold = (input.depth() == CV_32F ? FLT_EPSILON : DBL_EPSILON)*1000;
1621
1622     cvtest::convert( input, temp1, temp1.type() );
1623
1624     if( det < threshold ||
1625        ((method == CV_LU || method == CV_CHOLESKY) && (result == 0 || ratio < threshold)) ||
1626        ((method == CV_SVD || method == CV_SVD_SYM) && result < threshold) )
1627     {
1628         dst = Scalar::all(0);
1629         dst0 = Scalar::all(0);
1630         return;
1631     }
1632
1633     if( input.rows >= input.cols )
1634         cvtest::gemm( temp0, input, 1., Mat(), 0., dst, 0 );
1635     else
1636         cvtest::gemm( input, temp0, 1., Mat(), 0., dst, 0 );
1637
1638     cv::setIdentity( dst0, Scalar::all(1) );
1639 }
1640
1641
1642 ///////////////// solve /////////////////////
1643
1644 class Core_SolveTest : public Core_MatrixTest
1645 {
1646 public:
1647     typedef Core_MatrixTest Base;
1648     Core_SolveTest();
1649 protected:
1650     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1651     void get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high );
1652     double get_success_error_level( int test_case_idx, int i, int j );
1653     int prepare_test_case( int test_case_idx );
1654     void run_func();
1655     void prepare_to_validation( int test_case_idx );
1656     int method, rank;
1657     double result;
1658 };
1659
1660
1661 Core_SolveTest::Core_SolveTest() : Core_MatrixTest( 2, 1, false, false, 1 ), method(0), rank(0), result(0.)
1662 {
1663     test_case_count = 100;
1664     max_log_array_size = 7;
1665     test_array[TEMP].push_back(NULL);
1666     test_array[TEMP].push_back(NULL);
1667 }
1668
1669
1670 void Core_SolveTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
1671 {
1672     RNG& rng = ts->get_rng();
1673     int bits = cvtest::randInt(rng);
1674     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
1675     CvSize in_sz = sizes[INPUT][0];
1676     if( in_sz.width > in_sz.height )
1677         in_sz = cvSize(in_sz.height, in_sz.width);
1678     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
1679     sizes[INPUT][0] = in_sz;
1680     int min_size = MIN( sizes[INPUT][0].width, sizes[INPUT][0].height );
1681
1682     if( (bits & 3) == 0 )
1683     {
1684         method = CV_SVD;
1685         if( bits & 4 )
1686         {
1687             sizes[INPUT][0] = Size(min_size, min_size);
1688             /*if( bits & 8 )
1689              method = CV_SVD_SYM;*/
1690         }
1691     }
1692     else
1693     {
1694         method = CV_LU;
1695         sizes[INPUT][0] = Size(min_size, min_size);
1696     }
1697
1698     sizes[INPUT][1].height = sizes[INPUT][0].height;
1699     sizes[TEMP][0].width = sizes[INPUT][1].width;
1700     sizes[TEMP][0].height = sizes[INPUT][0].width;
1701     sizes[TEMP][1] = sizes[INPUT][0];
1702     types[TEMP][0] = types[INPUT][0];
1703     types[TEMP][1] = CV_64FC1;
1704     sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = Size(sizes[INPUT][1].width, min_size);
1705 }
1706
1707
1708 int Core_SolveTest::prepare_test_case( int test_case_idx )
1709 {
1710     int code = Core_MatrixTest::prepare_test_case( test_case_idx );
1711
1712     /*if( method == CV_SVD_SYM )
1713      {
1714      cvTsGEMM( test_array[INPUT][0], test_array[INPUT][0], 1.,
1715      0, 0., test_array[TEMP][0], CV_GEMM_B_T );
1716      cvTsCopy( test_array[TEMP][0], test_array[INPUT][0] );
1717      }*/
1718
1719     return code;
1720 }
1721
1722
1723 void Core_SolveTest::get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high )
1724 {
1725     low = cvScalarAll(-1.);
1726     high = cvScalarAll(1.);
1727 }
1728
1729
1730 double Core_SolveTest::get_success_error_level( int /*test_case_idx*/, int, int )
1731 {
1732     return CV_MAT_DEPTH(cvGetElemType(test_array[OUTPUT][0])) == CV_32F ? 5e-2 : 1e-8;
1733 }
1734
1735
1736 void Core_SolveTest::run_func()
1737 {
1738     result = cvSolve(test_array[INPUT][0], test_array[INPUT][1], test_array[TEMP][0], method);
1739 }
1740
1741 void Core_SolveTest::prepare_to_validation( int )
1742 {
1743     //int rank = test_mat[REF_OUTPUT][0].rows;
1744     Mat& input = test_mat[INPUT][0];
1745     Mat& dst = test_mat[OUTPUT][0];
1746     Mat& dst0 = test_mat[REF_OUTPUT][0];
1747
1748     if( method == CV_LU )
1749     {
1750         if( result == 0 )
1751         {
1752             Mat& temp1 = test_mat[TEMP][1];
1753             cvtest::convert(input, temp1, temp1.type());
1754             dst = Scalar::all(0);
1755             CvMat _temp1 = temp1;
1756             double det = cvTsLU( &_temp1, 0, 0 );
1757             dst0 = Scalar::all(det != 0);
1758             return;
1759         }
1760
1761         double threshold = (input.type() == CV_32F ? FLT_EPSILON : DBL_EPSILON)*1000;
1762         CvMat _input = input;
1763         double ratio = 0, det = cvTsSVDet( &_input, &ratio );
1764         if( det < threshold || ratio < threshold )
1765         {
1766             dst = Scalar::all(0);
1767             dst0 = Scalar::all(0);
1768             return;
1769         }
1770     }
1771
1772     Mat* pdst = input.rows <= input.cols ? &test_mat[OUTPUT][0] : &test_mat[INPUT][1];
1773
1774     cvtest::gemm( input, test_mat[TEMP][0], 1., test_mat[INPUT][1], -1., *pdst, 0 );
1775     if( pdst != &dst )
1776         cvtest::gemm( input, *pdst, 1., Mat(), 0., dst, CV_GEMM_A_T );
1777     dst0 = Scalar::all(0);
1778 }
1779
1780
1781 ///////////////// SVD /////////////////////
1782
1783 class Core_SVDTest : public Core_MatrixTest
1784 {
1785 public:
1786     typedef Core_MatrixTest Base;
1787     Core_SVDTest();
1788 protected:
1789     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
1790     double get_success_error_level( int test_case_idx, int i, int j );
1791     void get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high );
1792     int prepare_test_case( int test_case_idx );
1793     void run_func();
1794     void prepare_to_validation( int test_case_idx );
1795     int flags;
1796     bool have_u, have_v, symmetric, compact, vector_w;
1797 };
1798
1799
1800 Core_SVDTest::Core_SVDTest() :
1801 Core_MatrixTest( 1, 4, false, false, 1 ),
1802 flags(0), have_u(false), have_v(false), symmetric(false), compact(false), vector_w(false)
1803 {
1804     test_case_count = 100;
1805     test_array[TEMP].push_back(NULL);
1806     test_array[TEMP].push_back(NULL);
1807     test_array[TEMP].push_back(NULL);
1808     test_array[TEMP].push_back(NULL);
1809 }
1810
1811
1812 void Core_SVDTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types )
1813 {
1814     RNG& rng = ts->get_rng();
1815     int bits = cvtest::randInt(rng);
1816     Core_MatrixTest::get_test_array_types_and_sizes( test_case_idx, sizes, types );
1817     int min_size, i, m, n;
1818
1819     min_size = MIN( sizes[INPUT][0].width, sizes[INPUT][0].height );
1820
1821     flags = bits & (CV_SVD_MODIFY_A+CV_SVD_U_T+CV_SVD_V_T);
1822     have_u = (bits & 8) != 0;
1823     have_v = (bits & 16) != 0;
1824     symmetric = (bits & 32) != 0;
1825     compact = (bits & 64) != 0;
1826     vector_w = (bits & 128) != 0;
1827
1828     if( symmetric )
1829         sizes[INPUT][0] = Size(min_size, min_size);
1830
1831     m = sizes[INPUT][0].height;
1832     n = sizes[INPUT][0].width;
1833
1834     if( compact )
1835         sizes[TEMP][0] = Size(min_size, min_size);
1836     else
1837         sizes[TEMP][0] = sizes[INPUT][0];
1838     sizes[TEMP][3] = Size(0,0);
1839
1840     if( vector_w )
1841     {
1842         sizes[TEMP][3] = sizes[TEMP][0];
1843         if( bits & 256 )
1844             sizes[TEMP][0] = Size(1, min_size);
1845         else
1846             sizes[TEMP][0] = Size(min_size, 1);
1847     }
1848
1849     if( have_u )
1850     {
1851         sizes[TEMP][1] = compact ? Size(min_size, m) : Size(m, m);
1852
1853         if( flags & CV_SVD_U_T )
1854             CV_SWAP( sizes[TEMP][1].width, sizes[TEMP][1].height, i );
1855     }
1856     else
1857         sizes[TEMP][1] = Size(0,0);
1858
1859     if( have_v )
1860     {
1861         sizes[TEMP][2] = compact ? Size(n, min_size) : Size(n, n);
1862
1863         if( !(flags & CV_SVD_V_T) )
1864             CV_SWAP( sizes[TEMP][2].width, sizes[TEMP][2].height, i );
1865     }
1866     else
1867         sizes[TEMP][2] = Size(0,0);
1868
1869     types[TEMP][0] = types[TEMP][1] = types[TEMP][2] = types[TEMP][3] = types[INPUT][0];
1870     types[OUTPUT][0] = types[OUTPUT][1] = types[OUTPUT][2] = types[INPUT][0];
1871     types[OUTPUT][3] = CV_8UC1;
1872     sizes[OUTPUT][0] = !have_u || !have_v ? Size(0,0) : sizes[INPUT][0];
1873     sizes[OUTPUT][1] = !have_u ? Size(0,0) : compact ? Size(min_size,min_size) : Size(m,m);
1874     sizes[OUTPUT][2] = !have_v ? Size(0,0) : compact ? Size(min_size,min_size) : Size(n,n);
1875     sizes[OUTPUT][3] = Size(min_size,1);
1876
1877     for( i = 0; i < 4; i++ )
1878     {
1879         sizes[REF_OUTPUT][i] = sizes[OUTPUT][i];
1880         types[REF_OUTPUT][i] = types[OUTPUT][i];
1881     }
1882 }
1883
1884
1885 int Core_SVDTest::prepare_test_case( int test_case_idx )
1886 {
1887     int code = Core_MatrixTest::prepare_test_case( test_case_idx );
1888     if( code > 0 )
1889     {
1890         Mat& input = test_mat[INPUT][0];
1891         cvTsFloodWithZeros( input, ts->get_rng() );
1892
1893         if( symmetric && (have_u || have_v) )
1894         {
1895             Mat& temp = test_mat[TEMP][have_u ? 1 : 2];
1896             cvtest::gemm( input, input, 1., Mat(), 0., temp, CV_GEMM_B_T );
1897             cvtest::copy( temp, input );
1898         }
1899
1900         if( (flags & CV_SVD_MODIFY_A) && test_array[OUTPUT][0] )
1901             cvtest::copy( input, test_mat[OUTPUT][0] );
1902     }
1903
1904     return code;
1905 }
1906
1907
1908 void Core_SVDTest::get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high )
1909 {
1910     low = cvScalarAll(-2.);
1911     high = cvScalarAll(2.);
1912 }
1913
1914 double Core_SVDTest::get_success_error_level( int test_case_idx, int i, int j )
1915 {
1916     int input_depth = CV_MAT_DEPTH(cvGetElemType( test_array[INPUT][0] ));
1917     double input_precision = input_depth < CV_32F ? 0 : input_depth == CV_32F ? 1e-5 : 5e-11;
1918     double output_precision = Base::get_success_error_level( test_case_idx, i, j );
1919     return MAX(input_precision, output_precision);
1920 }
1921
1922 void Core_SVDTest::run_func()
1923 {
1924     CvArr* src = test_array[!(flags & CV_SVD_MODIFY_A) ? INPUT : OUTPUT][0];
1925     if( !src )
1926         src = test_array[INPUT][0];
1927     cvSVD( src, test_array[TEMP][0], test_array[TEMP][1], test_array[TEMP][2], flags );
1928 }
1929
1930
1931 void Core_SVDTest::prepare_to_validation( int /*test_case_idx*/ )
1932 {
1933     Mat& input = test_mat[INPUT][0];
1934     int depth = input.depth();
1935     int i, m = input.rows, n = input.cols, min_size = MIN(m, n);
1936     Mat *src, *dst, *w;
1937     double prev = 0, threshold = depth == CV_32F ? FLT_EPSILON : DBL_EPSILON;
1938
1939     if( have_u )
1940     {
1941         src = &test_mat[TEMP][1];
1942         dst = &test_mat[OUTPUT][1];
1943         cvtest::gemm( *src, *src, 1., Mat(), 0., *dst, src->rows == dst->rows ? CV_GEMM_B_T : CV_GEMM_A_T );
1944         cv::setIdentity( test_mat[REF_OUTPUT][1], Scalar::all(1.) );
1945     }
1946
1947     if( have_v )
1948     {
1949         src = &test_mat[TEMP][2];
1950         dst = &test_mat[OUTPUT][2];
1951         cvtest::gemm( *src, *src, 1., Mat(), 0., *dst, src->rows == dst->rows ? CV_GEMM_B_T : CV_GEMM_A_T );
1952         cv::setIdentity( test_mat[REF_OUTPUT][2], Scalar::all(1.) );
1953     }
1954
1955     w = &test_mat[TEMP][0];
1956     for( i = 0; i < min_size; i++ )
1957     {
1958         double normval = 0, aii;
1959         if( w->rows > 1 && w->cols > 1 )
1960         {
1961             normval = cvtest::norm( w->row(i), NORM_L1 );
1962             aii = depth == CV_32F ? w->at<float>(i,i) : w->at<double>(i,i);
1963         }
1964         else
1965         {
1966             normval = aii = depth == CV_32F ? w->at<float>(i) : w->at<double>(i);
1967         }
1968
1969         normval = fabs(normval - aii);
1970         test_mat[OUTPUT][3].at<uchar>(i) = aii >= 0 && normval < threshold && (i == 0 || aii <= prev);
1971         prev = aii;
1972     }
1973
1974     test_mat[REF_OUTPUT][3] = Scalar::all(1);
1975
1976     if( have_u && have_v )
1977     {
1978         if( vector_w )
1979         {
1980             test_mat[TEMP][3] = Scalar::all(0);
1981             for( i = 0; i < min_size; i++ )
1982             {
1983                 double val = depth == CV_32F ? w->at<float>(i) : w->at<double>(i);
1984                 cvSetReal2D( test_array[TEMP][3], i, i, val );
1985             }
1986             w = &test_mat[TEMP][3];
1987         }
1988
1989         if( m >= n )
1990         {
1991             cvtest::gemm( test_mat[TEMP][1], *w, 1., Mat(), 0., test_mat[REF_OUTPUT][0],
1992                      flags & CV_SVD_U_T ? CV_GEMM_A_T : 0 );
1993             cvtest::gemm( test_mat[REF_OUTPUT][0], test_mat[TEMP][2], 1., Mat(), 0.,
1994                      test_mat[OUTPUT][0], flags & CV_SVD_V_T ? 0 : CV_GEMM_B_T );
1995         }
1996         else
1997         {
1998             cvtest::gemm( *w, test_mat[TEMP][2], 1., Mat(), 0., test_mat[REF_OUTPUT][0],
1999                      flags & CV_SVD_V_T ? 0 : CV_GEMM_B_T );
2000             cvtest::gemm( test_mat[TEMP][1], test_mat[REF_OUTPUT][0], 1., Mat(), 0.,
2001                      test_mat[OUTPUT][0], flags & CV_SVD_U_T ? CV_GEMM_A_T : 0 );
2002         }
2003
2004         cvtest::copy( test_mat[INPUT][0], test_mat[REF_OUTPUT][0] );
2005     }
2006 }
2007
2008
2009
2010 ///////////////// SVBkSb /////////////////////
2011
2012 class Core_SVBkSbTest : public Core_MatrixTest
2013 {
2014 public:
2015     typedef Core_MatrixTest Base;
2016     Core_SVBkSbTest();
2017 protected:
2018     void get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes, vector<vector<int> >& types );
2019     double get_success_error_level( int test_case_idx, int i, int j );
2020     void get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high );
2021     int prepare_test_case( int test_case_idx );
2022     void run_func();
2023     void prepare_to_validation( int test_case_idx );
2024     int flags;
2025     bool have_b, symmetric, compact, vector_w;
2026 };
2027
2028
2029 Core_SVBkSbTest::Core_SVBkSbTest() : Core_MatrixTest( 2, 1, false, false, 1 ),
2030 flags(0), have_b(false), symmetric(false), compact(false), vector_w(false)
2031 {
2032     test_case_count = 100;
2033     test_array[TEMP].push_back(NULL);
2034     test_array[TEMP].push_back(NULL);
2035     test_array[TEMP].push_back(NULL);
2036 }
2037
2038
2039 void Core_SVBkSbTest::get_test_array_types_and_sizes( int test_case_idx, vector<vector<Size> >& sizes,
2040                                                       vector<vector<int> >& types )
2041 {
2042     RNG& rng = ts->get_rng();
2043     int bits = cvtest::randInt(rng);
2044     Base::get_test_array_types_and_sizes( test_case_idx, sizes, types );
2045     int min_size, i, m, n;
2046     CvSize b_size;
2047
2048     min_size = MIN( sizes[INPUT][0].width, sizes[INPUT][0].height );
2049
2050     flags = bits & (CV_SVD_MODIFY_A+CV_SVD_U_T+CV_SVD_V_T);
2051     have_b = (bits & 16) != 0;
2052     symmetric = (bits & 32) != 0;
2053     compact = (bits & 64) != 0;
2054     vector_w = (bits & 128) != 0;
2055
2056     if( symmetric )
2057         sizes[INPUT][0] = Size(min_size, min_size);
2058
2059     m = sizes[INPUT][0].height;
2060     n = sizes[INPUT][0].width;
2061
2062     sizes[INPUT][1] = Size(0,0);
2063     b_size = Size(m,m);
2064     if( have_b )
2065     {
2066         sizes[INPUT][1].height = sizes[INPUT][0].height;
2067         sizes[INPUT][1].width = cvtest::randInt(rng) % 100 + 1;
2068         b_size = sizes[INPUT][1];
2069     }
2070
2071     if( compact )
2072         sizes[TEMP][0] = Size(min_size, min_size);
2073     else
2074         sizes[TEMP][0] = sizes[INPUT][0];
2075
2076     if( vector_w )
2077     {
2078         if( bits & 256 )
2079             sizes[TEMP][0] = Size(1, min_size);
2080         else
2081             sizes[TEMP][0] = Size(min_size, 1);
2082     }
2083
2084     sizes[TEMP][1] = compact ? Size(min_size, m) : Size(m, m);
2085
2086     if( flags & CV_SVD_U_T )
2087         CV_SWAP( sizes[TEMP][1].width, sizes[TEMP][1].height, i );
2088
2089     sizes[TEMP][2] = compact ? Size(n, min_size) : Size(n, n);
2090
2091     if( !(flags & CV_SVD_V_T) )
2092         CV_SWAP( sizes[TEMP][2].width, sizes[TEMP][2].height, i );
2093
2094     types[TEMP][0] = types[TEMP][1] = types[TEMP][2] = types[INPUT][0];
2095     types[OUTPUT][0] = types[REF_OUTPUT][0] = types[INPUT][0];
2096     sizes[OUTPUT][0] = sizes[REF_OUTPUT][0] = Size( b_size.width, n );
2097 }
2098
2099
2100 int Core_SVBkSbTest::prepare_test_case( int test_case_idx )
2101 {
2102     int code = Base::prepare_test_case( test_case_idx );
2103     if( code > 0 )
2104     {
2105         Mat& input = test_mat[INPUT][0];
2106         cvTsFloodWithZeros( input, ts->get_rng() );
2107
2108         if( symmetric )
2109         {
2110             Mat& temp = test_mat[TEMP][1];
2111             cvtest::gemm( input, input, 1., Mat(), 0., temp, CV_GEMM_B_T );
2112             cvtest::copy( temp, input );
2113         }
2114
2115         CvMat _input = input;
2116         cvSVD( &_input, test_array[TEMP][0], test_array[TEMP][1], test_array[TEMP][2], flags );
2117     }
2118
2119     return code;
2120 }
2121
2122
2123 void Core_SVBkSbTest::get_minmax_bounds( int /*i*/, int /*j*/, int /*type*/, Scalar& low, Scalar& high )
2124 {
2125     low = cvScalarAll(-2.);
2126     high = cvScalarAll(2.);
2127 }
2128
2129
2130 double Core_SVBkSbTest::get_success_error_level( int /*test_case_idx*/, int /*i*/, int /*j*/ )
2131 {
2132     return CV_MAT_DEPTH(cvGetElemType(test_array[INPUT][0])) == CV_32F ? 1e-3 : 1e-7;
2133 }
2134
2135
2136 void Core_SVBkSbTest::run_func()
2137 {
2138     cvSVBkSb( test_array[TEMP][0], test_array[TEMP][1], test_array[TEMP][2],
2139              test_array[INPUT][1], test_array[OUTPUT][0], flags );
2140 }
2141
2142
2143 void Core_SVBkSbTest::prepare_to_validation( int )
2144 {
2145     Mat& input = test_mat[INPUT][0];
2146     int i, m = input.rows, n = input.cols, min_size = MIN(m, n);
2147     bool is_float = input.type() == CV_32F;
2148     Size w_size = compact ? Size(min_size,min_size) : Size(m,n);
2149     Mat& w = test_mat[TEMP][0];
2150     Mat wdb( w_size.height, w_size.width, CV_64FC1 );
2151     CvMat _w = w, _wdb = wdb;
2152     // use exactly the same threshold as in icvSVD... ,
2153     // so the changes in the library and here should be synchronized.
2154     double threshold = cv::sum(w)[0]*(DBL_EPSILON*2);//(is_float ? FLT_EPSILON*10 : DBL_EPSILON*2);
2155
2156     wdb = Scalar::all(0);
2157     for( i = 0; i < min_size; i++ )
2158     {
2159         double wii = vector_w ? cvGetReal1D(&_w,i) : cvGetReal2D(&_w,i,i);
2160         cvSetReal2D( &_wdb, i, i, wii > threshold ? 1./wii : 0. );
2161     }
2162
2163     Mat u = test_mat[TEMP][1];
2164     Mat v = test_mat[TEMP][2];
2165     Mat b = test_mat[INPUT][1];
2166
2167     if( is_float )
2168     {
2169         test_mat[TEMP][1].convertTo(u, CV_64F);
2170         test_mat[TEMP][2].convertTo(v, CV_64F);
2171         if( !b.empty() )
2172             test_mat[INPUT][1].convertTo(b, CV_64F);
2173     }
2174
2175     Mat t0, t1;
2176
2177     if( !b.empty() )
2178         cvtest::gemm( u, b, 1., Mat(), 0., t0, !(flags & CV_SVD_U_T) ? CV_GEMM_A_T : 0 );
2179     else if( flags & CV_SVD_U_T )
2180         cvtest::copy( u, t0 );
2181     else
2182         cvtest::transpose( u, t0 );
2183
2184     cvtest::gemm( wdb, t0, 1, Mat(), 0, t1, 0 );
2185
2186     cvtest::gemm( v, t1, 1, Mat(), 0, t0, flags & CV_SVD_V_T ? CV_GEMM_A_T : 0 );
2187     Mat& dst0 = test_mat[REF_OUTPUT][0];
2188     t0.convertTo(dst0, dst0.type() );
2189 }
2190
2191
2192 typedef std::complex<double> complex_type;
2193
2194 struct pred_complex
2195 {
2196     bool operator() (const complex_type& lhs, const complex_type& rhs) const
2197     {
2198         return fabs(lhs.real() - rhs.real()) > fabs(rhs.real())*FLT_EPSILON ? lhs.real() < rhs.real() : lhs.imag() < rhs.imag();
2199     }
2200 };
2201
2202 struct pred_double
2203 {
2204     bool operator() (const double& lhs, const double& rhs) const
2205     {
2206         return lhs < rhs;
2207     }
2208 };
2209
2210 class Core_SolvePolyTest : public cvtest::BaseTest
2211 {
2212 public:
2213     Core_SolvePolyTest();
2214     ~Core_SolvePolyTest();
2215 protected:
2216     virtual void run( int start_from );
2217 };
2218
2219 Core_SolvePolyTest::Core_SolvePolyTest() {}
2220
2221 Core_SolvePolyTest::~Core_SolvePolyTest() {}
2222
2223 void Core_SolvePolyTest::run( int )
2224 {
2225     RNG& rng = ts->get_rng();
2226     int fig = 100;
2227     double range = 50;
2228     double err_eps = 1e-4;
2229
2230     for (int idx = 0, max_idx = 1000, progress = 0; idx < max_idx; ++idx)
2231     {
2232         progress = update_progress(progress, idx-1, max_idx, 0);
2233         int n = cvtest::randInt(rng) % 13 + 1;
2234         std::vector<complex_type> r(n), ar(n), c(n + 1, 0);
2235         std::vector<double> a(n + 1), u(n * 2), ar1(n), ar2(n);
2236
2237         int rr_odds = 3; // odds that we get a real root
2238         for (int j = 0; j < n;)
2239         {
2240             if (cvtest::randInt(rng) % rr_odds == 0 || j == n - 1)
2241                     r[j++] = cvtest::randReal(rng) * range;
2242             else
2243             {
2244                     r[j] = complex_type(cvtest::randReal(rng) * range,
2245                                     cvtest::randReal(rng) * range + 1);
2246                     r[j + 1] = std::conj(r[j]);
2247                     j += 2;
2248             }
2249         }
2250
2251         for (int j = 0, k = 1 << n, jj, kk; j < k; ++j)
2252         {
2253             int p = 0;
2254             complex_type v(1);
2255             for (jj = 0, kk = 1; jj < n && !(j & kk); ++jj, ++p, kk <<= 1)
2256                 ;
2257             for (; jj < n; ++jj, kk <<= 1)
2258             {
2259                     if (j & kk)
2260                         v *= -r[jj];
2261                     else
2262                         ++p;
2263             }
2264             c[p] += v;
2265         }
2266
2267         bool pass = false;
2268         double div = 0, s = 0;
2269         int cubic_case = idx & 1;
2270         for (int maxiter = 100; !pass && maxiter < 10000; maxiter *= 2, cubic_case = (cubic_case + 1) % 2)
2271         {
2272             for (int j = 0; j < n + 1; ++j)
2273                     a[j] = c[j].real();
2274
2275             CvMat amat, umat;
2276             cvInitMatHeader(&amat, n + 1, 1, CV_64FC1, &a[0]);
2277             cvInitMatHeader(&umat, n, 1, CV_64FC2, &u[0]);
2278             cvSolvePoly(&amat, &umat, maxiter, fig);
2279
2280             for (int j = 0; j < n; ++j)
2281                     ar[j] = complex_type(u[j * 2], u[j * 2 + 1]);
2282
2283             std::sort(r.begin(), r.end(), pred_complex());
2284             std::sort(ar.begin(), ar.end(), pred_complex());
2285
2286             pass = true;
2287             if( n == 3 )
2288             {
2289                 ar2.resize(n);
2290                 cv::Mat _umat2(3, 1, CV_64F, &ar2[0]), umat2 = _umat2;
2291                 cvFlip(&amat, &amat, 0);
2292                 int nr2;
2293                 if( cubic_case == 0 )
2294                     nr2 = cv::solveCubic(cv::Mat(&amat),umat2);
2295                 else
2296                     nr2 = cv::solveCubic(cv::Mat_<float>(cv::Mat(&amat)), umat2);
2297                 cvFlip(&amat, &amat, 0);
2298                 if(nr2 > 0)
2299                     std::sort(ar2.begin(), ar2.begin()+nr2, pred_double());
2300                 ar2.resize(nr2);
2301
2302                 int nr1 = 0;
2303                 for(int j = 0; j < n; j++)
2304                     if( fabs(r[j].imag()) < DBL_EPSILON )
2305                         ar1[nr1++] = r[j].real();
2306
2307                 pass = pass && nr1 == nr2;
2308                 if( nr2 > 0 )
2309                 {
2310                     div = s = 0;
2311                     for(int j = 0; j < nr1; j++)
2312                     {
2313                         s += fabs(ar1[j]);
2314                         div += fabs(ar1[j] - ar2[j]);
2315                     }
2316                     div /= s;
2317                     pass = pass && div < err_eps;
2318                 }
2319             }
2320
2321             div = s = 0;
2322             for (int j = 0; j < n; ++j)
2323             {
2324                 s += fabs(r[j].real()) + fabs(r[j].imag());
2325                 div += sqrt(pow(r[j].real() - ar[j].real(), 2) + pow(r[j].imag() - ar[j].imag(), 2));
2326             }
2327             div /= s;
2328             pass = pass && div < err_eps;
2329         }
2330
2331         if (!pass)
2332         {
2333             ts->set_failed_test_info(cvtest::TS::FAIL_INVALID_OUTPUT);
2334             ts->printf( cvtest::TS::LOG, "too big diff = %g\n", div );
2335
2336             for (size_t j=0;j<ar2.size();++j)
2337                 ts->printf( cvtest::TS::LOG, "ar2[%d]=%g\n", j, ar2[j]);
2338             ts->printf(cvtest::TS::LOG, "\n");
2339
2340             for (size_t j=0;j<r.size();++j)
2341                     ts->printf( cvtest::TS::LOG, "r[%d]=(%g, %g)\n", j, r[j].real(), r[j].imag());
2342             ts->printf( cvtest::TS::LOG, "\n" );
2343             for (size_t j=0;j<ar.size();++j)
2344                     ts->printf( cvtest::TS::LOG, "ar[%d]=(%g, %g)\n", j, ar[j].real(), ar[j].imag());
2345             break;
2346         }
2347     }
2348 }
2349
2350
2351 /////////////////////////////////////////////////////////////////////////////////////////////////////
2352
2353 TEST(Core_CovarMatrix, accuracy) { Core_CovarMatrixTest test; test.safe_run(); }
2354 TEST(Core_CrossProduct, accuracy) { Core_CrossProductTest test; test.safe_run(); }
2355 TEST(Core_Determinant, accuracy) { Core_DetTest test; test.safe_run(); }
2356 TEST(Core_DotProduct, accuracy) { Core_DotProductTest test; test.safe_run(); }
2357 TEST(Core_GEMM, accuracy) { Core_GEMMTest test; test.safe_run(); }
2358 TEST(Core_Invert, accuracy) { Core_InvertTest test; test.safe_run(); }
2359 TEST(Core_Mahalanobis, accuracy) { Core_MahalanobisTest test; test.safe_run(); }
2360 TEST(Core_MulTransposed, accuracy) { Core_MulTransposedTest test; test.safe_run(); }
2361 TEST(Core_Transform, accuracy) { Core_TransformTest test; test.safe_run(); }
2362 TEST(Core_PerspectiveTransform, accuracy) { Core_PerspectiveTransformTest test; test.safe_run(); }
2363 TEST(Core_Pow, accuracy) { Core_PowTest test; test.safe_run(); }
2364 TEST(Core_SolveLinearSystem, accuracy) { Core_SolveTest test; test.safe_run(); }
2365 TEST(Core_SVD, accuracy) { Core_SVDTest test; test.safe_run(); }
2366 TEST(Core_SVBkSb, accuracy) { Core_SVBkSbTest test; test.safe_run(); }
2367 TEST(Core_Trace, accuracy) { Core_TraceTest test; test.safe_run(); }
2368 TEST(Core_SolvePoly, accuracy) { Core_SolvePolyTest test; test.safe_run(); }
2369
2370 // TODO: eigenvv, invsqrt, cbrt, fastarctan, (round, floor, ceil(?)),
2371
2372 /* End of file. */
2373