modules/core/test/test_rand.cpp

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