• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
Core_RandTest()17 Core_RandTest::Core_RandTest()
18 {
19 }
20 
chi2_p95(int n)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 
check_pdf(const Mat & hist,double scale,int dist_type,double & refval,double & realval)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 
run(int)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 
TEST(Core_Rand,quality)305 TEST(Core_Rand, quality) { Core_RandTest test; test.safe_run(); }
306 
307 
308 class Core_RandRangeTest : public cvtest::BaseTest
309 {
310 public:
Core_RandRangeTest()311     Core_RandRangeTest() {}
~Core_RandRangeTest()312     ~Core_RandRangeTest() {}
313 protected:
run(int)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 
TEST(Core_Rand,range)340 TEST(Core_Rand, range) { Core_RandRangeTest test; test.safe_run(); }
341 
342 
TEST(Core_RNG_MT19937,regression)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