// This file is part of OpenCV project. // It is subject to the license terms in the LICENSE file found in the top-level directory // of this distribution and at http://opencv.org/license.html. #include "test_precomp.hpp" namespace opencv_test { namespace { class Core_RandTest : public cvtest::BaseTest { public: Core_RandTest(); protected: void run(int); bool check_pdf(const Mat& hist, double scale, int dist_type, double& refval, double& realval); }; Core_RandTest::Core_RandTest() { } static double chi2_p95(int n) { static float chi2_tab95[] = { 3.841f, 5.991f, 7.815f, 9.488f, 11.07f, 12.59f, 14.07f, 15.51f, 16.92f, 18.31f, 19.68f, 21.03f, 21.03f, 22.36f, 23.69f, 25.00f, 26.30f, 27.59f, 28.87f, 30.14f, 31.41f, 32.67f, 33.92f, 35.17f, 36.42f, 37.65f, 38.89f, 40.11f, 41.34f, 42.56f, 43.77f }; static const double xp = 1.64; CV_Assert(n >= 1); if( n <= 30 ) return chi2_tab95[n-1]; return n + sqrt((double)2*n)*xp + 0.6666666666666*(xp*xp - 1); } bool Core_RandTest::check_pdf(const Mat& hist, double scale, int dist_type, double& refval, double& realval) { Mat hist0(hist.size(), CV_32F); const int* H = hist.ptr<int>(); float* H0 = hist0.ptr<float>(); int i, hsz = hist.cols; double sum = 0; for( i = 0; i < hsz; i++ ) sum += H[i]; CV_Assert( fabs(1./sum - scale) < FLT_EPSILON ); if( dist_type == CV_RAND_UNI ) { float scale0 = (float)(1./hsz); for( i = 0; i < hsz; i++ ) H0[i] = scale0; } else { double sum2 = 0, r = (hsz-1.)/2; double alpha = 2*sqrt(2.)/r, beta = -alpha*r; for( i = 0; i < hsz; i++ ) { double x = i*alpha + beta; H0[i] = (float)exp(-x*x); sum2 += H0[i]; } sum2 = 1./sum2; for( i = 0; i < hsz; i++ ) H0[i] = (float)(H0[i]*sum2); } double chi2 = 0; for( i = 0; i < hsz; i++ ) { double a = H0[i]; double b = H[i]*scale; if( a > DBL_EPSILON ) chi2 += (a - b)*(a - b)/(a + b); } realval = chi2; double chi2_pval = chi2_p95(hsz - 1 - (dist_type == CV_RAND_NORMAL ? 2 : 0)); refval = chi2_pval*0.01; return realval <= refval; } void Core_RandTest::run( int ) { static int _ranges[][2] = {{ 0, 256 }, { -128, 128 }, { 0, 65536 }, { -32768, 32768 }, { -1000000, 1000000 }, { -1000, 1000 }, { -1000, 1000 }}; const int MAX_SDIM = 10; const int N = 2000000; const int maxSlice = 1000; const int MAX_HIST_SIZE = 1000; int progress = 0; RNG& rng = ts->get_rng(); RNG tested_rng = theRNG(); test_case_count = 200; for( int idx = 0; idx < test_case_count; idx++ ) { progress = update_progress( progress, idx, test_case_count, 0 ); ts->update_context( this, idx, false ); int depth = cvtest::randInt(rng) % (CV_64F+1); int c, cn = (cvtest::randInt(rng) % 4) + 1; int type = CV_MAKETYPE(depth, cn); int dist_type = cvtest::randInt(rng) % (CV_RAND_NORMAL+1); int i, k, SZ = N/cn; Scalar A, B; double eps = 1.e-4; if (depth == CV_64F) eps = 1.e-7; bool do_sphere_test = dist_type == CV_RAND_UNI; Mat arr[2], hist[4]; int W[] = {0,0,0,0}; arr[0].create(1, SZ, type); arr[1].create(1, SZ, type); bool fast_algo = dist_type == CV_RAND_UNI && depth < CV_32F; for( c = 0; c < cn; c++ ) { int a, b, hsz; if( dist_type == CV_RAND_UNI ) { a = (int)(cvtest::randInt(rng) % (_ranges[depth][1] - _ranges[depth][0])) + _ranges[depth][0]; do { b = (int)(cvtest::randInt(rng) % (_ranges[depth][1] - _ranges[depth][0])) + _ranges[depth][0]; } while( abs(a-b) <= 1 ); if( a > b ) std::swap(a, b); unsigned r = (unsigned)(b - a); fast_algo = fast_algo && r <= 256 && (r & (r-1)) == 0; hsz = min((unsigned)(b - a), (unsigned)MAX_HIST_SIZE); do_sphere_test = do_sphere_test && b - a >= 100; } else { int vrange = _ranges[depth][1] - _ranges[depth][0]; int meanrange = vrange/16; int mindiv = MAX(vrange/20, 5); int maxdiv = MIN(vrange/8, 10000); a = cvtest::randInt(rng) % meanrange - meanrange/2 + (_ranges[depth][0] + _ranges[depth][1])/2; b = cvtest::randInt(rng) % (maxdiv - mindiv) + mindiv; hsz = min((unsigned)b*9, (unsigned)MAX_HIST_SIZE); } A[c] = a; B[c] = b; hist[c].create(1, hsz, CV_32S); } cv::RNG saved_rng = tested_rng; int maxk = fast_algo ? 0 : 1; for( k = 0; k <= maxk; k++ ) { tested_rng = saved_rng; int sz = 0, dsz = 0, slice; for( slice = 0; slice < maxSlice && sz < SZ; slice++, sz += dsz ) { dsz = slice+1 < maxSlice ? (int)(cvtest::randInt(rng) % (SZ - sz) + 1) : SZ - sz; Mat aslice = arr[k].colRange(sz, sz + dsz); tested_rng.fill(aslice, dist_type, A, B); } } if( maxk >= 1 && cvtest::norm(arr[0], arr[1], NORM_INF) > eps) { ts->printf( cvtest::TS::LOG, "RNG output depends on the array lengths (some generated numbers get lost?)" ); ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT ); return; } for( c = 0; c < cn; c++ ) { const uchar* data = arr[0].ptr(); int* H = hist[c].ptr<int>(); int HSZ = hist[c].cols; double minVal = dist_type == CV_RAND_UNI ? A[c] : A[c] - B[c]*4; double maxVal = dist_type == CV_RAND_UNI ? B[c] : A[c] + B[c]*4; double scale = HSZ/(maxVal - minVal); double delta = -minVal*scale; hist[c] = Scalar::all(0); for( i = c; i < SZ*cn; i += cn ) { double val = depth == CV_8U ? ((const uchar*)data)[i] : depth == CV_8S ? ((const schar*)data)[i] : depth == CV_16U ? ((const ushort*)data)[i] : depth == CV_16S ? ((const short*)data)[i] : depth == CV_32S ? ((const int*)data)[i] : depth == CV_32F ? ((const float*)data)[i] : ((const double*)data)[i]; int ival = cvFloor(val*scale + delta); if( (unsigned)ival < (unsigned)HSZ ) { H[ival]++; W[c]++; } else if( dist_type == CV_RAND_UNI ) { if( (minVal <= val && val < maxVal) || (depth >= CV_32F && val == maxVal) ) { H[ival < 0 ? 0 : HSZ-1]++; W[c]++; } else { putchar('^'); } } } if( dist_type == CV_RAND_UNI && W[c] != SZ ) { ts->printf( cvtest::TS::LOG, "Uniform RNG gave values out of the range [%g,%g) on channel %d/%d\n", A[c], B[c], c, cn); ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT ); return; } if( dist_type == CV_RAND_NORMAL && W[c] < SZ*.90) { 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", A[c], B[c], A[c], B[c], c, cn); ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT ); return; } double refval = 0, realval = 0; if( !check_pdf(hist[c], 1./W[c], dist_type, refval, realval) ) { ts->printf( cvtest::TS::LOG, "RNG failed Chi-square test " "(got %g vs probable maximum %g) on channel %d/%d\n", realval, refval, c, cn); ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT ); return; } } // Monte-Carlo test. Compute volume of SDIM-dimensional sphere // inscribed in [-1,1]^SDIM cube. if( do_sphere_test ) { int SDIM = cvtest::randInt(rng) % (MAX_SDIM-1) + 2; int N0 = (SZ*cn/SDIM), n = 0; double r2 = 0; const uchar* data = arr[0].ptr(); double scale[4], delta[4]; for( c = 0; c < cn; c++ ) { scale[c] = 2./(B[c] - A[c]); delta[c] = -A[c]*scale[c] - 1; } for( i = k = c = 0; i <= SZ*cn - SDIM; i++, k++, c++ ) { double val = depth == CV_8U ? ((const uchar*)data)[i] : depth == CV_8S ? ((const schar*)data)[i] : depth == CV_16U ? ((const ushort*)data)[i] : depth == CV_16S ? ((const short*)data)[i] : depth == CV_32S ? ((const int*)data)[i] : depth == CV_32F ? ((const float*)data)[i] : ((const double*)data)[i]; c &= c < cn ? -1 : 0; val = val*scale[c] + delta[c]; r2 += val*val; if( k == SDIM-1 ) { n += r2 <= 1; r2 = 0; k = -1; } } double V = ((double)n/N0)*(1 << SDIM); // the theoretically computed volume int sdim = SDIM % 2; double V0 = sdim + 1; for( sdim += 2; sdim <= SDIM; sdim += 2 ) V0 *= 2*CV_PI/sdim; if( fabs(V - V0) > 0.3*fabs(V0) ) { ts->printf( cvtest::TS::LOG, "RNG failed %d-dim sphere volume test (got %g instead of %g)\n", SDIM, V, V0); ts->printf( cvtest::TS::LOG, "depth = %d, N0 = %d\n", depth, N0); ts->set_failed_test_info( cvtest::TS::FAIL_INVALID_OUTPUT ); return; } } } } TEST(Core_Rand, quality) { Core_RandTest test; test.safe_run(); } class Core_RandRangeTest : public cvtest::BaseTest { public: Core_RandRangeTest() {} ~Core_RandRangeTest() {} protected: void run(int) { Mat a(Size(1280, 720), CV_8U, Scalar(20)); Mat af(Size(1280, 720), CV_32F, Scalar(20)); theRNG().fill(a, RNG::UNIFORM, -DBL_MAX, DBL_MAX); theRNG().fill(af, RNG::UNIFORM, -DBL_MAX, DBL_MAX); int n0 = 0, n255 = 0, nx = 0; int nfmin = 0, nfmax = 0, nfx = 0; for( int i = 0; i < a.rows; i++ ) for( int j = 0; j < a.cols; j++ ) { int v = a.at<uchar>(i,j); double vf = af.at<float>(i,j); if( v == 0 ) n0++; else if( v == 255 ) n255++; else nx++; if( vf < FLT_MAX*-0.999f ) nfmin++; else if( vf > FLT_MAX*0.999f ) nfmax++; else nfx++; } CV_Assert( n0 > nx*2 && n255 > nx*2 ); CV_Assert( nfmin > nfx*2 && nfmax > nfx*2 ); } }; TEST(Core_Rand, range) { Core_RandRangeTest test; test.safe_run(); } TEST(Core_RNG_MT19937, regression) { cv::RNG_MT19937 rng; int actual[61] = {0, }; const size_t length = (sizeof(actual) / sizeof(actual[0])); for (int i = 0; i < 10000; ++i ) { actual[(unsigned)(rng.next() ^ i) % length]++; } int expected[length] = { 177, 158, 180, 177, 160, 179, 143, 162, 177, 144, 170, 174, 165, 168, 168, 156, 177, 157, 159, 169, 177, 182, 166, 154, 144, 180, 168, 152, 170, 187, 160, 145, 139, 164, 157, 179, 148, 183, 159, 160, 196, 184, 149, 142, 162, 148, 163, 152, 168, 173, 160, 181, 172, 181, 155, 153, 158, 171, 138, 150, 150 }; for (size_t i = 0; i < length; ++i) { ASSERT_EQ(expected[i], actual[i]); } } TEST(Core_Rand, Regression_Stack_Corruption) { int bufsz = 128; //enough for 14 doubles AutoBuffer<uchar> buffer(bufsz); size_t offset = 0; cv::Mat_<cv::Point2d> x(2, 3, (cv::Point2d*)(buffer.data()+offset)); offset += x.total()*x.elemSize(); double& param1 = *(double*)(buffer.data()+offset); offset += sizeof(double); double& param2 = *(double*)(buffer.data()+offset); param1 = -9; param2 = 2; cv::theRNG().fill(x, cv::RNG::NORMAL, param1, param2); ASSERT_EQ(param1, -9); ASSERT_EQ(param2, 2); } class RandRowFillParallelLoopBody : public cv::ParallelLoopBody { public: RandRowFillParallelLoopBody(Mat& dst) : dst_(dst) {} ~RandRowFillParallelLoopBody() {} void operator()(const cv::Range& r) const { cv::RNG rng = cv::theRNG(); // copy state for (int y = r.start; y < r.end; y++) { cv::theRNG() = cv::RNG(rng.state + y); // seed is based on processed row cv::randu(dst_.row(y), Scalar(-100), Scalar(100)); } // theRNG() state is changed here (but state collision has low probability, so we don't check this) } protected: Mat& dst_; }; TEST(Core_Rand, parallel_for_stable_results) { cv::RNG rng = cv::theRNG(); // save rng state Mat dst1(1000, 100, CV_8SC1); parallel_for_(cv::Range(0, dst1.rows), RandRowFillParallelLoopBody(dst1)); cv::theRNG() = rng; // restore rng state Mat dst2(1000, 100, CV_8SC1); parallel_for_(cv::Range(0, dst2.rows), RandRowFillParallelLoopBody(dst2)); ASSERT_EQ(0, countNonZero(dst1 != dst2)); } }} // namespace