opencv/tests/cxcore/src/arand.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41
  42 #include "cxcoretest.h"
  43
  44 using namespace cv;
  45
  46
  47 class CV_RandTest : public CvTest
  48 {
  49 public:
  50     CV_RandTest();
  51 protected:
  52     void run(int);
  53     bool check_pdf(const Mat& hist, double scale, double A, double B,
  54                    int dist_type, double& refval, double& realval);
  55 };
  56
  57
  58 CV_RandTest::CV_RandTest():
  59 CvTest( "rand", "cvRandArr, cvRNG" )
  60 {
  61     support_testing_modes = CvTS::CORRECTNESS_CHECK_MODE;
  62 }
  63
  64 static double chi2_p95(int n)
  65 {
  66     static float chi2_tab95[] = {
  67         3.841, 5.991, 7.815, 9.488, 11.07, 12.59, 14.07, 15.51, 16.92, 18.31, 19.68, 21.03,
  68         21.03, 22.36, 23.69, 25.00, 26.30, 27.59, 28.87, 30.14, 31.41, 32.67, 33.92, 35.17,
  69         36.42, 37.65, 38.89, 40.11, 41.34, 42.56, 43.77 };
  70     static const double xp = 1.64;
  71     CV_Assert(n >= 1);
  72
  73     if( n <= 30 )
  74         return chi2_tab95[n-1];
  75     return n + sqrt((double)2*n)*xp + 0.6666666666666*(xp*xp - 1);
  76 }
  77
  78 bool CV_RandTest::check_pdf(const Mat& hist, double scale, double A, double B,
  79                             int dist_type, double& refval, double& realval)
  80 {
  81     Mat hist0(hist.size(), CV_32F);
  82     const int* H = (const int*)hist.data;
  83     float* H0 = ((float*)hist0.data);
  84     int i, hsz = hist.cols;
  85
  86     double sum = 0;
  87     for( i = 0; i < hsz; i++ )
  88         sum += H[i];
  89     CV_Assert( fabs(1./sum - scale) < FLT_EPSILON );
  90
  91     if( dist_type == CV_RAND_UNI )
  92     {
  93         float scale0 = (float)(1./hsz);
  94         for( i = 0; i < hsz; i++ )
  95             H0[i] = scale0;
  96     }
  97     else
  98     {
  99         double sum = 0, r = (hsz-1.)/2;
 100         double alpha = 2*sqrt(2.)/r, beta = -alpha*r;
 101         for( i = 0; i < hsz; i++ )
 102         {
 103             double x = i*alpha + beta;
 104             H0[i] = (float)exp(-x*x);
 105             sum += H0[i];
 106         }
 107         sum = 1./sum;
 108         for( i = 0; i < hsz; i++ )
 109             H0[i] = (float)(H0[i]*sum);
 110     }
 111
 112     double chi2 = 0;
 113     for( i = 0; i < hsz; i++ )
 114     {
 115         double a = H0[i];
 116         double b = H[i]*scale;
 117         if( a > DBL_EPSILON )
 118             chi2 += (a - b)*(a - b)/(a + b);
 119     }
 120
 121     double chi2_pval = chi2_p95(hsz - 1 - (dist_type == CV_RAND_NORMAL ? 2 : 0));
 122     return chi2 <= chi2_pval*0.01;
 123 }
 124
 125 void CV_RandTest::run( int start_from )
 126 {
 127     static int _ranges[][2] =
 128     {{ 0, 256 }, { -128, 128 }, { 0, 65536 }, { -32768, 32768 },
 129         { -1000000, 1000000 }, { -1000, 1000 }, { -1000, 1000 }};
 130
 131     const int MAX_SDIM = 10;
 132     const int N = 1200000;
 133     const int maxSlice = 1000;
 134     const int MAX_HIST_SIZE = 1000;
 135     int progress = 0;
 136
 137     CvRNG* rng = ts->get_rng();
 138     RNG tested_rng;
 139     test_case_count = 500;
 140
 141     for( int idx = 0; idx < test_case_count; idx++ )
 142     {
 143         ts->update_context( this, idx, false );
 144
 145         int depth = cvTsRandInt(rng) % (CV_64F+1);
 146         int c, cn = (cvTsRandInt(rng) % 4) + 1;
 147         int type = CV_MAKETYPE(depth, cn);
 148         int dist_type = cvTsRandInt(rng) % (CV_RAND_NORMAL+1);
 149         int i, k, SZ = N/cn;
 150         Scalar A, B;
 151
 152         bool do_sphere_test = dist_type == CV_RAND_UNI;
 153         Mat arr[2], hist[4];
 154         int W[] = {0,0,0,0};
 155
 156         arr[0].create(1, SZ, type);
 157         arr[1].create(1, SZ, type);
 158         bool fast_algo = dist_type == CV_RAND_UNI && depth < CV_32F;
 159
 160         for( c = 0; c < cn; c++ )
 161         {
 162             int a, b, hsz;
 163             if( dist_type == CV_RAND_UNI )
 164             {
 165                 a = (int)(cvTsRandInt(rng) % (_ranges[depth][1] -
 166                         _ranges[depth][0])) + _ranges[depth][0];
 167                 do
 168                 {
 169                     b = (int)(cvTsRandInt(rng) % (_ranges[depth][1] -
 170                         _ranges[depth][0])) + _ranges[depth][0];
 171                 }
 172                 while( abs(a-b) <= 1 );
 173                 if( a > b )
 174                     std::swap(a, b);
 175
 176                 unsigned r = (unsigned)(b - a);
 177                 fast_algo = fast_algo && r <= 256 && (r & (r-1)) == 0;
 178                 hsz = min((unsigned)(b - a), (unsigned)MAX_HIST_SIZE);
 179                 do_sphere_test = do_sphere_test && b - a >= 100;
 180             }
 181             else
 182             {
 183                 int vrange = _ranges[depth][1] - _ranges[depth][0];
 184                 int meanrange = vrange/16;
 185                 int mindiv = MAX(vrange/20, 5);
 186                 int maxdiv = MIN(vrange/8, 10000);
 187
 188                 a = cvTsRandInt(rng) % meanrange - meanrange/2 +
 189                               (_ranges[depth][0] + _ranges[depth][1])/2;
 190                 b = cvTsRandInt(rng) % (maxdiv - mindiv) + mindiv;
 191                 hsz = min((unsigned)b*9, (unsigned)MAX_HIST_SIZE);
 192             }
 193             A[c] = a;
 194             B[c] = b;
 195             hist[c].create(1, hsz, CV_32S);
 196         }
 197
 198         cv::RNG saved_rng = tested_rng;
 199         int maxk = fast_algo ? 0 : 1;
 200         for( k = 0; k <= maxk; k++ )
 201         {
 202             tested_rng = saved_rng;
 203             int sz = 0, dsz = 0, slice;
 204             for( slice = 0; slice < maxSlice; slice++, sz += dsz )
 205             {
 206                 dsz = slice+1 < maxSlice ? cvTsRandInt(rng) % (SZ - sz + 1) : SZ - sz;
 207                 Mat aslice = arr[k].colRange(sz, sz + dsz);
 208                 tested_rng.fill(aslice, dist_type, A, B);
 209             }
 210         }
 211
 212         if( maxk >= 1 && norm(arr[0], arr[1], NORM_INF) != 0 )
 213         {
 214             ts->printf( CvTS::LOG, "RNG output depends on the array lengths (some generated numbers get lost?)" );
 215             ts->set_failed_test_info( CvTS::FAIL_INVALID_OUTPUT );
 216             return;
 217         }
 218
 219         for( c = 0; c < cn; c++ )
 220         {
 221             const uchar* data = arr[0].data;
 222             int* H = hist[c].ptr<int>();
 223             int HSZ = hist[c].cols;
 224             double minVal = dist_type == CV_RAND_UNI ? A[c] : A[c] - B[c]*4;
 225             double maxVal = dist_type == CV_RAND_UNI ? B[c] : A[c] + B[c]*4;
 226             double scale = HSZ/(maxVal - minVal);
 227             double delta = -minVal*scale;
 228
 229             hist[c] = Scalar::all(0);
 230
 231             for( i = c; i < SZ*cn; i += cn )
 232             {
 233                 double val = depth == CV_8U ? ((const uchar*)data)[i] :
 234                     depth == CV_8S ? ((const schar*)data)[i] :
 235                     depth == CV_16U ? ((const ushort*)data)[i] :
 236                     depth == CV_16S ? ((const short*)data)[i] :
 237                     depth == CV_32S ? ((const int*)data)[i] :
 238                     depth == CV_32F ? ((const float*)data)[i] :
 239                                       ((const double*)data)[i];
 240                 int ival = cvFloor(val*scale + delta);
 241                 if( (unsigned)ival < (unsigned)HSZ )
 242                 {
 243                     H[ival]++;
 244                     W[c]++;
 245                 }
 246                 else if( dist_type == CV_RAND_UNI )
 247                 {
 248                     if( depth >= CV_32F && val == maxVal )
 249                     {
 250                         H[HSZ-1]++;
 251                         W[c]++;
 252                     }
 253                     else
 254                     {
 255                         putchar('^');
 256                     }
 257                 }
 258             }
 259
 260             if( dist_type == CV_RAND_UNI && W[c] != SZ )
 261             {
 262                 ts->printf( CvTS::LOG, "Uniform RNG gave values out of the range [%g,%g) on channel %d/%d\n",
 263                            A[c], B[c], c, cn);
 264                 ts->set_failed_test_info( CvTS::FAIL_INVALID_OUTPUT );
 265                 return;
 266             }
 267             if( dist_type == CV_RAND_NORMAL && W[c] < SZ*.90)
 268             {
 269                 ts->printf( CvTS::LOG, "Normal RNG gave too many values out of the range (%g+4*%g,%g+4*%g) on channel %d/%d\n",
 270                            A[c], B[c], A[c], B[c], c, cn);
 271                 ts->set_failed_test_info( CvTS::FAIL_INVALID_OUTPUT );
 272                 return;
 273             }
 274             double refval = 0, realval = 0;
 275
 276             if( !check_pdf(hist[c], 1./W[c], A[c], B[c], dist_type, refval, realval) )
 277             {
 278                 ts->printf( CvTS::LOG, "RNG failed Chi-square test "
 279                            "(got %g vs probable maximum %g) on channel %d/%d\n",
 280                            realval, refval, c, cn);
 281                 ts->set_failed_test_info( CvTS::FAIL_INVALID_OUTPUT );
 282                 return;
 283             }
 284         }
 285
 286         // Monte-Carlo test. Compute volume of SDIM-dimensional sphere
 287         // inscribed in [-1,1]^SDIM cube.
 288         if( do_sphere_test )
 289         {
 290             int SDIM = cvTsRandInt(rng) % (MAX_SDIM-1) + 2;
 291             int N0 = (SZ*cn/SDIM), N = 0;
 292             double r2 = 0;
 293             const uchar* data = arr[0].data;
 294             double scale[4], delta[4];
 295             for( c = 0; c < cn; c++ )
 296             {
 297                 scale[c] = 2./(B[c] - A[c]);
 298                 delta[c] = -A[c]*scale[c] - 1;
 299             }
 300
 301             for( i = k = c = 0; i <= SZ*cn - SDIM; i++, k++, c++ )
 302             {
 303                 double val = depth == CV_8U ? ((const uchar*)data)[i] :
 304                     depth == CV_8S ? ((const schar*)data)[i] :
 305                     depth == CV_16U ? ((const ushort*)data)[i] :
 306                     depth == CV_16S ? ((const short*)data)[i] :
 307                     depth == CV_32S ? ((const int*)data)[i] :
 308                     depth == CV_32F ? ((const float*)data)[i] : ((const double*)data)[i];
 309                 c &= c < cn ? -1 : 0;
 310                 val = val*scale[c] + delta[c];
 311                 r2 += val*val;
 312                 if( k == SDIM-1 )
 313                 {
 314                     N += r2 <= 1;
 315                     r2 = 0;
 316                     k = -1;
 317                 }
 318             }
 319
 320             double V = ((double)N/N0)*(1 << SDIM);
 321
 322             // the theoretically computed volume
 323             int sdim = SDIM % 2;
 324             double V0 = sdim + 1;
 325             for( sdim += 2; sdim <= SDIM; sdim += 2 )
 326                 V0 *= 2*CV_PI/sdim;
 327
 328             if( fabs(V - V0) > 0.2*fabs(V0) )
 329             {
 330                 ts->printf( CvTS::LOG, "RNG failed %d-dim sphere volume test (got %g instead of %g)\n",
 331                            SDIM, V, V0);
 332                 ts->set_failed_test_info( CvTS::FAIL_INVALID_OUTPUT );
 333                 return;
 334             }
 335         }
 336         progress = update_progress( progress, idx, test_case_count, 0 );
 337     }
 338 }
 339
 340 CV_RandTest rand_test;