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 = (const int*)hist.data;
  41     float* H0 = ((float*)hist0.data);
  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 ? 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 && 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].data;
 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].data;
 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