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
