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 && sz < SZ; 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                 printf("%d - %d\n", sz, sz + dsz);
 177             }
 178         }
 179
 180         if( maxk >= 1 && cvtest::norm(arr[0], arr[1], NORM_INF) > eps)
 181         {
 182             ts->printf( cvtest::TS::LOG, "RNG output depends on the array lengths (some generated numbers get lost?)" );
 183             ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 184             return;
 185         }
 186
 187         for( c = 0; c < cn; c++ )
 188         {
 189             const uchar* data = arr[0].ptr();
 190             int* H = hist[c].ptr<int>();
 191             int HSZ = hist[c].cols;
 192             double minVal = dist_type == CV_RAND_UNI ? A[c] : A[c] - B[c]*4;
 193             double maxVal = dist_type == CV_RAND_UNI ? B[c] : A[c] + B[c]*4;
 194             double scale = HSZ/(maxVal - minVal);
 195             double delta = -minVal*scale;
 196
 197             hist[c] = Scalar::all(0);
 198
 199             for( i = c; i < SZ*cn; i += cn )
 200             {
 201                 double val = depth == CV_8U ? ((const uchar*)data)[i] :
 202                 depth == CV_8S ? ((const schar*)data)[i] :
 203                 depth == CV_16U ? ((const ushort*)data)[i] :
 204                 depth == CV_16S ? ((const short*)data)[i] :
 205                 depth == CV_32S ? ((const int*)data)[i] :
 206                 depth == CV_32F ? ((const float*)data)[i] :
 207                 ((const double*)data)[i];
 208                 int ival = cvFloor(val*scale + delta);
 209                 if( (unsigned)ival < (unsigned)HSZ )
 210                 {
 211                     H[ival]++;
 212                     W[c]++;
 213                 }
 214                 else if( dist_type == CV_RAND_UNI )
 215                 {
 216                     if( (minVal <= val && val < maxVal) || (depth >= CV_32F && val == maxVal) )
 217                     {
 218                         H[ival < 0 ? 0 : HSZ-1]++;
 219                         W[c]++;
 220                     }
 221                     else
 222                     {
 223                         putchar('^');
 224                     }
 225                 }
 226             }
 227
 228             if( dist_type == CV_RAND_UNI && W[c] != SZ )
 229             {
 230                 ts->printf( cvtest::TS::LOG, "Uniform RNG gave values out of the range [%g,%g) on channel %d/%d\n",
 231                            A[c], B[c], c, cn);
 232                 ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 233                 return;
 234             }
 235             if( dist_type == CV_RAND_NORMAL && W[c] < SZ*.90)
 236             {
 237                 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",
 238                            A[c], B[c], A[c], B[c], c, cn);
 239                 ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 240                 return;
 241             }
 242             double refval = 0, realval = 0;
 243
 244             if( !check_pdf(hist[c], 1./W[c], dist_type, refval, realval) )
 245             {
 246                 ts->printf( cvtest::TS::LOG, "RNG failed Chi-square test "
 247                            "(got %g vs probable maximum %g) on channel %d/%d\n",
 248                            realval, refval, c, cn);
 249                 ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 250                 return;
 251             }
 252         }
 253
 254         // Monte-Carlo test. Compute volume of SDIM-dimensional sphere
 255         // inscribed in [-1,1]^SDIM cube.
 256         if( do_sphere_test )
 257         {
 258             int SDIM = cvtest::randInt(rng) % (MAX_SDIM-1) + 2;
 259             int N0 = (SZ*cn/SDIM), n = 0;
 260             double r2 = 0;
 261             const uchar* data = arr[0].ptr();
 262             double scale[4], delta[4];
 263             for( c = 0; c < cn; c++ )
 264             {
 265                 scale[c] = 2./(B[c] - A[c]);
 266                 delta[c] = -A[c]*scale[c] - 1;
 267             }
 268
 269             for( i = k = c = 0; i <= SZ*cn - SDIM; i++, k++, c++ )
 270             {
 271                 double val = depth == CV_8U ? ((const uchar*)data)[i] :
 272                 depth == CV_8S ? ((const schar*)data)[i] :
 273                 depth == CV_16U ? ((const ushort*)data)[i] :
 274                 depth == CV_16S ? ((const short*)data)[i] :
 275                 depth == CV_32S ? ((const int*)data)[i] :
 276                 depth == CV_32F ? ((const float*)data)[i] : ((const double*)data)[i];
 277                 c &= c < cn ? -1 : 0;
 278                 val = val*scale[c] + delta[c];
 279                 r2 += val*val;
 280                 if( k == SDIM-1 )
 281                 {
 282                     n += r2 <= 1;
 283                     r2 = 0;
 284                     k = -1;
 285                 }
 286             }
 287
 288             double V = ((double)n/N0)*(1 << SDIM);
 289
 290             // the theoretically computed volume
 291             int sdim = SDIM % 2;
 292             double V0 = sdim + 1;
 293             for( sdim += 2; sdim <= SDIM; sdim += 2 )
 294                 V0 *= 2*CV_PI/sdim;
 295
 296             if( fabs(V - V0) > 0.3*fabs(V0) )
 297             {
 298                 ts->printf( cvtest::TS::LOG, "RNG failed %d-dim sphere volume test (got %g instead of %g)\n",
 299                            SDIM, V, V0);
 300                 ts->printf( cvtest::TS::LOG, "depth = %d, N0 = %d\n", depth, N0);
 301                 ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT );
 302                 return;
 303             }
 304         }
 305     }
 306 }
 307
 308 TEST(Core_Rand, quality) { Core_RandTest test; test.safe_run(); }
 309
 310
 311 class Core_RandRangeTest : public cvtest::BaseTest
 312 {
 313 public:
 314     Core_RandRangeTest() {}
 315     ~Core_RandRangeTest() {}
 316 protected:
 317     void run(int)
 318     {
 319         Mat a(Size(1280, 720), CV_8U, Scalar(20));
 320         Mat af(Size(1280, 720), CV_32F, Scalar(20));
 321         theRNG().fill(a, RNG::UNIFORM, -DBL_MAX, DBL_MAX);
 322         theRNG().fill(af, RNG::UNIFORM, -DBL_MAX, DBL_MAX);
 323         int n0 = 0, n255 = 0, nx = 0;
 324         int nfmin = 0, nfmax = 0, nfx = 0;
 325
 326         for( int i = 0; i < a.rows; i++ )
 327             for( int j = 0; j < a.cols; j++ )
 328             {
 329                 int v = a.at<uchar>(i,j);
 330                 double vf = af.at<float>(i,j);
 331                 if( v == 0 ) n0++;
 332                 else if( v == 255 ) n255++;
 333                 else nx++;
 334                 if( vf < FLT_MAX*-0.999f ) nfmin++;
 335                 else if( vf > FLT_MAX*0.999f ) nfmax++;
 336                 else nfx++;
 337             }
 338         CV_Assert( n0 > nx*2 && n255 > nx*2 );
 339         CV_Assert( nfmin > nfx*2 && nfmax > nfx*2 );
 340     }
 341 };
 342
 343 TEST(Core_Rand, range) { Core_RandRangeTest test; test.safe_run(); }
 344
 345
 346 TEST(Core_RNG_MT19937, regression)
 347 {
 348     cv::RNG_MT19937 rng;
 349     int actual[61] = {0, };
 350     const size_t length = (sizeof(actual) / sizeof(actual[0]));
 351     for (int i = 0; i < 10000; ++i )
 352     {
 353         actual[(unsigned)(rng.next() ^ i) % length]++;
 354     }
 355
 356     int expected[length] = {
 357         177, 158, 180, 177,  160, 179, 143, 162,
 358         177, 144, 170, 174,  165, 168, 168, 156,
 359         177, 157, 159, 169,  177, 182, 166, 154,
 360         144, 180, 168, 152,  170, 187, 160, 145,
 361         139, 164, 157, 179,  148, 183, 159, 160,
 362         196, 184, 149, 142,  162, 148, 163, 152,
 363         168, 173, 160, 181,  172, 181, 155, 153,
 364         158, 171, 138, 150,  150 };
 365
 366     for (size_t i = 0; i < length; ++i)
 367     {
 368         ASSERT_EQ(expected[i], actual[i]);
 369     }
 370 }
 371
 372
 373 TEST(Core_Rand, Regression_Stack_Corruption)
 374 {
 375     int bufsz = 128; //enough for 14 doubles
 376     AutoBuffer<uchar> buffer(bufsz);
 377     size_t offset = 0;
 378     cv::Mat_<cv::Point2d> x(2, 3, (cv::Point2d*)(buffer.data()+offset)); offset += x.total()*x.elemSize();
 379     double& param1 = *(double*)(buffer.data()+offset); offset += sizeof(double);
 380     double& param2 = *(double*)(buffer.data()+offset); offset += sizeof(double);
 381     param1 = -9; param2 = 2;
 382
 383     cv::theRNG().fill(x, cv::RNG::NORMAL, param1, param2);
 384
 385     ASSERT_EQ(param1, -9);
 386     ASSERT_EQ(param2,  2);
 387 }
 388
 389
 390 class RandRowFillParallelLoopBody : public cv::ParallelLoopBody
 391 {
 392 public:
 393     RandRowFillParallelLoopBody(Mat& dst) : dst_(dst) {}
 394     ~RandRowFillParallelLoopBody() {}
 395     void operator()(const cv::Range& r) const
 396     {
 397         cv::RNG rng = cv::theRNG(); // copy state
 398         for (int y = r.start; y < r.end; y++)
 399         {
 400             cv::theRNG() = cv::RNG(rng.state + y); // seed is based on processed row
 401             cv::randu(dst_.row(y), Scalar(-100), Scalar(100));
 402         }
 403         // theRNG() state is changed here (but state collision has low probability, so we don't check this)
 404     }
 405 protected:
 406     Mat& dst_;
 407 };
 408
 409 TEST(Core_Rand, parallel_for_stable_results)
 410 {
 411     cv::RNG rng = cv::theRNG(); // save rng state
 412     Mat dst1(1000, 100, CV_8SC1);
 413     parallel_for_(cv::Range(0, dst1.rows), RandRowFillParallelLoopBody(dst1));
 414
 415     cv::theRNG() = rng; // restore rng state
 416     Mat dst2(1000, 100, CV_8SC1);
 417     parallel_for_(cv::Range(0, dst2.rows), RandRowFillParallelLoopBody(dst2));
 418
 419     ASSERT_EQ(0, countNonZero(dst1 != dst2));
 420 }
 421
 422 }} // namespace