modules/core/test/test_rand.cpp

   1 // This file is part of OpenCV project.
   2 // It is subject to the license terms in the LICENSE file found in the top-level directory
   3 // of this distribution and at http://opencv.org/license.html.
   4 #include "test_precomp.hpp"
   5
   6 namespace opencv_test { namespace {
   7
   8 class Core_RandTest : public cvtest::BaseTest
   9 {
  10 public:
  11     Core_RandTest();
  12 protected:
  13     void run(int);
  14     bool check_pdf(const Mat& hist, double scale, int dist_type,
  15                    double& refval, double& realval);
  16 };
  17
  18
  19 Core_RandTest::Core_RandTest()
  20 {
  21 }
  22
  23 static double chi2_p95(int n)
  24 {
  25     static float chi2_tab95[] = {
  26         3.841f, 5.991f, 7.815f, 9.488f, 11.07f, 12.59f, 14.07f, 15.51f,
  27         16.92f, 18.31f, 19.68f, 21.03f, 21.03f, 22.36f, 23.69f, 25.00f,
  28         26.30f, 27.59f, 28.87f, 30.14f, 31.41f, 32.67f, 33.92f, 35.17f,
  29         36.42f, 37.65f, 38.89f, 40.11f, 41.34f, 42.56f, 43.77f };
  30     static const double xp = 1.64;
  31     CV_Assert(n >= 1);
  32
  33     if( n <= 30 )
  34         return chi2_tab95[n-1];
  35     return n + sqrt((double)2*n)*xp + 0.6666666666666*(xp*xp - 1);
  36 }
  37
  38 bool Core_RandTest::check_pdf(const Mat& hist, double scale,
  39                             int dist_type, double& refval, double& realval)
  40 {
  41     Mat hist0(hist.size(), CV_32F);
  42     const int* H = hist.ptr<int>();
  43     float* H0 = hist0.ptr<float>();
  44     int i, hsz = hist.cols;
  45
  46     double sum = 0;
  47     for( i = 0; i < hsz; i++ )
  48         sum += H[i];
  49     CV_Assert( fabs(1./sum - scale) < FLT_EPSILON );
  50
  51     if( dist_type == CV_RAND_UNI )
  52     {
  53         float scale0 = (float)(1./hsz);
  54         for( i = 0; i < hsz; i++ )
  55             H0[i] = scale0;
  56     }
  57     else
  58     {
  59         double sum2 = 0, r = (hsz-1.)/2;
  60         double alpha = 2*sqrt(2.)/r, beta = -alpha*r;
  61         for( i = 0; i < hsz; i++ )
  62         {
  63             double x = i*alpha + beta;
  64             H0[i] = (float)exp(-x*x);
  65             sum2 += H0[i];
  66         }
  67         sum2 = 1./sum2;
  68         for( i = 0; i < hsz; i++ )
  69             H0[i] = (float)(H0[i]*sum2);
  70     }
  71
  72     double chi2 = 0;
  73     for( i = 0; i < hsz; i++ )
  74     {
  75         double a = H0[i];
  76         double b = H[i]*scale;
  77         if( a > DBL_EPSILON )
  78             chi2 += (a - b)*(a - b)/(a + b);
  79     }
  80     realval = chi2;
  81
  82     double chi2_pval = chi2_p95(hsz - 1 - (dist_type == CV_RAND_NORMAL ? 2 : 0));
  83     refval = chi2_pval*0.01;
  84     return realval <= refval;
  85 }
  86
  87 void Core_RandTest::run( int )
  88 {
  89     static int _ranges[][2] =
  90     {{ 0, 256 }, { -128, 128 }, { 0, 65536 }, { -32768, 32768 },
  91         { -1000000, 1000000 }, { -1000, 1000 }, { -1000, 1000 }};
  92
  93     const int MAX_SDIM = 10;
  94     const int N = 2000000;
  95     const int maxSlice = 1000;
  96     const int MAX_HIST_SIZE = 1000;
  97     int progress = 0;
  98
  99     RNG& rng = ts->get_rng();
 100     RNG tested_rng = theRNG();
 101     test_case_count = 200;
 102
 103     for( int idx = 0; idx < test_case_count; idx++ )
 104     {
 105         progress = update_progress( progress, idx, test_case_count, 0 );
 106         ts->update_context( this, idx, false );
 107
 108         int depth = cvtest::randInt(rng) % (CV_64F+1);
 109         int c, cn = (cvtest::randInt(rng) % 4) + 1;
 110         int type = CV_MAKETYPE(depth, cn);
 111         int dist_type = cvtest::randInt(rng) % (CV_RAND_NORMAL+1);
 112         int i, k, SZ = N/cn;
 113         Scalar A, B;
 114
 115         double eps = 1.e-4;
 116         if (depth == CV_64F)
 117             eps = 1.e-7;
 118
 119         bool do_sphere_test = dist_type == CV_RAND_UNI;
 120         Mat arr[2], hist[4];
 121         int W[] = {0,0,0,0};
 122
 123         arr[0].create(1, SZ, type);
 124         arr[1].create(1, SZ, type);
 125         bool fast_algo = dist_type == CV_RAND_UNI && depth < CV_32F;
 126
 127         for( c = 0; c < cn; c++ )
 128         {
 129             int a, b, hsz;
 130             if( dist_type == CV_RAND_UNI )
 131             {
 132                 a = (int)(cvtest::randInt(rng) % (_ranges[depth][1] -
 133                                               _ranges[depth][0])) + _ranges[depth][0];
 134                 do
 135                 {
 136                     b = (int)(cvtest::randInt(rng) % (_ranges[depth][1] -
 137                                                   _ranges[depth][0])) + _ranges[depth][0];
 138                 }
 139                 while( abs(a-b) <= 1 );
 140                 if( a > b )
 141                     std::swap(a, b);
 142
 143                 unsigned r = (unsigned)(b - a);
 144                 fast_algo = fast_algo && r <= 256 && (r & (r-1)) == 0;
 145                 hsz = min((unsigned)(b - a), (unsigned)MAX_HIST_SIZE);
 146                 do_sphere_test = do_sphere_test && b - a >= 100;
 147             }
 148             else
 149             {
 150                 int vrange = _ranges[depth][1] - _ranges[depth][0];
 151                 int meanrange = vrange/16;
 152                 int mindiv = MAX(vrange/20, 5);
 153                 int maxdiv = MIN(vrange/8, 10000);
 154
 155                 a = cvtest::randInt(rng) % meanrange - meanrange/2 +
 156                 (_ranges[depth][0] + _ranges[depth][1])/2;
 157                 b = cvtest::randInt(rng) % (maxdiv - mindiv) + mindiv;
 158                 hsz = min((unsigned)b*9, (unsigned)MAX_HIST_SIZE);
 159             }
 160             A[c] = a;
 161             B[c] = b;
 162             hist[c].create(1, hsz, CV_32S);
 163         }
 164
 165         cv::RNG saved_rng = tested_rng;
 166         int maxk = fast_algo ? 0 : 1;
 167         for( k = 0; k <= maxk; k++ )
 168         {
 169             tested_rng = saved_rng;
 170             int sz = 0, dsz = 0, slice;
 171             for( slice = 0; slice < maxSlice; slice++, sz += dsz )
 172             {
 173                 dsz = slice+1 < maxSlice ? (int)(cvtest::randInt(rng) % (SZ - sz + 1)) : SZ - sz;
 174                 Mat aslice = arr[k].colRange(sz, sz + dsz);
 175                 tested_rng.fill(aslice, dist_type, A, B);
 176             }
 177         }
 178
 179         if( maxk >= 1 && cvtest::norm(arr[0], arr[1], NORM_INF) > eps)
 180         {
 181             ts->printf( cvtest::TS::LOG, "RNG output depends on the array lengths (some generated numbers get lost?)" );
 182             ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 183             return;
 184         }
 185
 186         for( c = 0; c < cn; c++ )
 187         {
 188             const uchar* data = arr[0].ptr();
 189             int* H = hist[c].ptr<int>();
 190             int HSZ = hist[c].cols;
 191             double minVal = dist_type == CV_RAND_UNI ? A[c] : A[c] - B[c]*4;
 192             double maxVal = dist_type == CV_RAND_UNI ? B[c] : A[c] + B[c]*4;
 193             double scale = HSZ/(maxVal - minVal);
 194             double delta = -minVal*scale;
 195
 196             hist[c] = Scalar::all(0);
 197
 198             for( i = c; i < SZ*cn; i += cn )
 199             {
 200                 double val = depth == CV_8U ? ((const uchar*)data)[i] :
 201                 depth == CV_8S ? ((const schar*)data)[i] :
 202                 depth == CV_16U ? ((const ushort*)data)[i] :
 203                 depth == CV_16S ? ((const short*)data)[i] :
 204                 depth == CV_32S ? ((const int*)data)[i] :
 205                 depth == CV_32F ? ((const float*)data)[i] :
 206                 ((const double*)data)[i];
 207                 int ival = cvFloor(val*scale + delta);
 208                 if( (unsigned)ival < (unsigned)HSZ )
 209                 {
 210                     H[ival]++;
 211                     W[c]++;
 212                 }
 213                 else if( dist_type == CV_RAND_UNI )
 214                 {
 215                     if( (minVal <= val && val < maxVal) || (depth >= CV_32F && val == maxVal) )
 216                     {
 217                         H[ival < 0 ? 0 : HSZ-1]++;
 218                         W[c]++;
 219                     }
 220                     else
 221                     {
 222                         putchar('^');
 223                     }
 224                 }
 225             }
 226
 227             if( dist_type == CV_RAND_UNI && W[c] != SZ )
 228             {
 229                 ts->printf( cvtest::TS::LOG, "Uniform RNG gave values out of the range [%g,%g) on channel %d/%d\n",
 230                            A[c], B[c], c, cn);
 231                 ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 232                 return;
 233             }
 234             if( dist_type == CV_RAND_NORMAL && W[c] < SZ*.90)
 235             {
 236                 ts->printf( cvtest::TS::LOG, "Normal RNG gave too many values out of the range (%g+4*%g,%g+4*%g) on channel %d/%d\n",
 237                            A[c], B[c], A[c], B[c], c, cn);
 238                 ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 239                 return;
 240             }
 241             double refval = 0, realval = 0;
 242
 243             if( !check_pdf(hist[c], 1./W[c], dist_type, refval, realval) )
 244             {
 245                 ts->printf( cvtest::TS::LOG, "RNG failed Chi-square test "
 246                            "(got %g vs probable maximum %g) on channel %d/%d\n",
 247                            realval, refval, c, cn);
 248                 ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 249                 return;
 250             }
 251         }
 252
 253         // Monte-Carlo test. Compute volume of SDIM-dimensional sphere
 254         // inscribed in [-1,1]^SDIM cube.
 255         if( do_sphere_test )
 256         {
 257             int SDIM = cvtest::randInt(rng) % (MAX_SDIM-1) + 2;
 258             int N0 = (SZ*cn/SDIM), n = 0;
 259             double r2 = 0;
 260             const uchar* data = arr[0].ptr();
 261             double scale[4], delta[4];
 262             for( c = 0; c < cn; c++ )
 263             {
 264                 scale[c] = 2./(B[c] - A[c]);
 265                 delta[c] = -A[c]*scale[c] - 1;
 266             }
 267
 268             for( i = k = c = 0; i <= SZ*cn - SDIM; i++, k++, c++ )
 269             {
 270                 double val = depth == CV_8U ? ((const uchar*)data)[i] :
 271                 depth == CV_8S ? ((const schar*)data)[i] :
 272                 depth == CV_16U ? ((const ushort*)data)[i] :
 273                 depth == CV_16S ? ((const short*)data)[i] :
 274                 depth == CV_32S ? ((const int*)data)[i] :
 275                 depth == CV_32F ? ((const float*)data)[i] : ((const double*)data)[i];
 276                 c &= c < cn ? -1 : 0;
 277                 val = val*scale[c] + delta[c];
 278                 r2 += val*val;
 279                 if( k == SDIM-1 )
 280                 {
 281                     n += r2 <= 1;
 282                     r2 = 0;
 283                     k = -1;
 284                 }
 285             }
 286
 287             double V = ((double)n/N0)*(1 << SDIM);
 288
 289             // the theoretically computed volume
 290             int sdim = SDIM % 2;
 291             double V0 = sdim + 1;
 292             for( sdim += 2; sdim <= SDIM; sdim += 2 )
 293                 V0 *= 2*CV_PI/sdim;
 294
 295             if( fabs(V - V0) > 0.3*fabs(V0) )
 296             {
 297                 ts->printf( cvtest::TS::LOG, "RNG failed %d-dim sphere volume test (got %g instead of %g)\n",
 298                            SDIM, V, V0);
 299                 ts->printf( cvtest::TS::LOG, "depth = %d, N0 = %d\n", depth, N0);
 300                 ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 301                 return;
 302             }
 303         }
 304     }
 305 }
 306
 307 TEST(Core_Rand, quality) { Core_RandTest test; test.safe_run(); }
 308
 309
 310 class Core_RandRangeTest : public cvtest::BaseTest
 311 {
 312 public:
 313     Core_RandRangeTest() {}
 314     ~Core_RandRangeTest() {}
 315 protected:
 316     void run(int)
 317     {
 318         Mat a(Size(1280, 720), CV_8U, Scalar(20));
 319         Mat af(Size(1280, 720), CV_32F, Scalar(20));
 320         theRNG().fill(a, RNG::UNIFORM, -DBL_MAX, DBL_MAX);
 321         theRNG().fill(af, RNG::UNIFORM, -DBL_MAX, DBL_MAX);
 322         int n0 = 0, n255 = 0, nx = 0;
 323         int nfmin = 0, nfmax = 0, nfx = 0;
 324
 325         for( int i = 0; i < a.rows; i++ )
 326             for( int j = 0; j < a.cols; j++ )
 327             {
 328                 int v = a.at<uchar>(i,j);
 329                 double vf = af.at<float>(i,j);
 330                 if( v == 0 ) n0++;
 331                 else if( v == 255 ) n255++;
 332                 else nx++;
 333                 if( vf < FLT_MAX*-0.999f ) nfmin++;
 334                 else if( vf > FLT_MAX*0.999f ) nfmax++;
 335                 else nfx++;
 336             }
 337         CV_Assert( n0 > nx*2 && n255 > nx*2 );
 338         CV_Assert( nfmin > nfx*2 && nfmax > nfx*2 );
 339     }
 340 };
 341
 342 TEST(Core_Rand, range) { Core_RandRangeTest test; test.safe_run(); }
 343
 344
 345 TEST(Core_RNG_MT19937, regression)
 346 {
 347     cv::RNG_MT19937 rng;
 348     int actual[61] = {0, };
 349     const size_t length = (sizeof(actual) / sizeof(actual[0]));
 350     for (int i = 0; i < 10000; ++i )
 351     {
 352         actual[(unsigned)(rng.next() ^ i) % length]++;
 353     }
 354
 355     int expected[length] = {
 356         177, 158, 180, 177,  160, 179, 143, 162,
 357         177, 144, 170, 174,  165, 168, 168, 156,
 358         177, 157, 159, 169,  177, 182, 166, 154,
 359         144, 180, 168, 152,  170, 187, 160, 145,
 360         139, 164, 157, 179,  148, 183, 159, 160,
 361         196, 184, 149, 142,  162, 148, 163, 152,
 362         168, 173, 160, 181,  172, 181, 155, 153,
 363         158, 171, 138, 150,  150 };
 364
 365     for (size_t i = 0; i < length; ++i)
 366     {
 367         ASSERT_EQ(expected[i], actual[i]);
 368     }
 369 }
 370
 371
 372 TEST(Core_Rand, Regression_Stack_Corruption)
 373 {
 374     int bufsz = 128; //enough for 14 doubles
 375     AutoBuffer<uchar> buffer(bufsz);
 376     size_t offset = 0;
 377     cv::Mat_<cv::Point2d> x(2, 3, (cv::Point2d*)(buffer+offset)); offset += x.total()*x.elemSize();
 378     double& param1 = *(double*)(buffer+offset); offset += sizeof(double);
 379     double& param2 = *(double*)(buffer+offset); offset += sizeof(double);
 380     param1 = -9; param2 = 2;
 381
 382     cv::theRNG().fill(x, cv::RNG::NORMAL, param1, param2);
 383
 384     ASSERT_EQ(param1, -9);
 385     ASSERT_EQ(param2,  2);
 386 }
 387
 388
 389 class RandRowFillParallelLoopBody : public cv::ParallelLoopBody
 390 {
 391 public:
 392     RandRowFillParallelLoopBody(Mat& dst) : dst_(dst) {}
 393     ~RandRowFillParallelLoopBody() {}
 394     void operator()(const cv::Range& r) const
 395     {
 396         cv::RNG rng = cv::theRNG(); // copy state
 397         for (int y = r.start; y < r.end; y++)
 398         {
 399             cv::theRNG() = cv::RNG(rng.state + y); // seed is based on processed row
 400             cv::randu(dst_.row(y), Scalar(-100), Scalar(100));
 401         }
 402         // theRNG() state is changed here (but state collision has low probability, so we don't check this)
 403     }
 404 protected:
 405     Mat& dst_;
 406 };
 407
 408 TEST(Core_Rand, parallel_for_stable_results)
 409 {
 410     cv::RNG rng = cv::theRNG(); // save rng state
 411     Mat dst1(1000, 100, CV_8SC1);
 412     parallel_for_(cv::Range(0, dst1.rows), RandRowFillParallelLoopBody(dst1));
 413
 414     cv::theRNG() = rng; // restore rng state
 415     Mat dst2(1000, 100, CV_8SC1);
 416     parallel_for_(cv::Range(0, dst2.rows), RandRowFillParallelLoopBody(dst2));
 417
 418     ASSERT_EQ(0, countNonZero(dst1 != dst2));
 419 }
 420
 421 }} // namespace