--- /dev/null
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+// By downloading, copying, installing or using the software you agree to this license.
+// If you do not agree to this license, do not download, install,
+// copy or use the software.
+//
+//
+// License Agreement
+// For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+// * Redistribution's of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+//
+// * Redistribution's in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+//
+// * The name of the copyright holders may not be used to endorse or promote products
+// derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+/* ////////////////////////////////////////////////////////////////////
+//
+// Filling CvMat/IplImage instances with random numbers
+//
+// */
+
+#include "_cxcore.h"
+
+namespace cv
+{
+
+///////////////////////////// Functions Declaration //////////////////////////////////////
+
+/*
+ Multiply-with-carry generator is used here:
+ temp = ( A*X(n) + carry )
+ X(n+1) = temp mod (2^32)
+ carry = temp / (2^32)
+*/
+
+#define RNG_NEXT(x) ((uint64)(unsigned)(x)*RNG::A + ((x) >> 32))
+
+/***************************************************************************************\
+* Pseudo-Random Number Generators (PRNGs) *
+\***************************************************************************************/
+
+template<typename T> static void
+RandBits_( Mat& _arr, uint64* state, const void* _param )
+{
+ uint64 temp = *state;
+ const int* param = (const int*)_param;
+ int small_flag = (param[12]|param[13]|param[14]|param[15]) <= 255;
+ Size size = getContinuousSize(_arr,_arr.channels());
+
+ for( int y = 0; y < size.height; y++ )
+ {
+ T* arr = (T*)(_arr.data + _arr.step*y);
+ int i, k = 3;
+ const int* p = param;
+
+ if( !small_flag )
+ {
+ for( i = 0; i <= size.width - 4; i += 4 )
+ {
+ int t0, t1;
+
+ temp = RNG_NEXT(temp);
+ t0 = ((int)temp & p[i + 12]) + p[i];
+ temp = RNG_NEXT(temp);
+ t1 = ((int)temp & p[i + 13]) + p[i+1];
+ arr[i] = saturate_cast<T>(t0);
+ arr[i+1] = saturate_cast<T>(t1);
+
+ temp = RNG_NEXT(temp);
+ t0 = ((int)temp & p[i + 14]) + p[i+2];
+ temp = RNG_NEXT(temp);
+ t1 = ((int)temp & p[i + 15]) + p[i+3];
+ arr[i+2] = saturate_cast<T>(t0);
+ arr[i+3] = saturate_cast<T>(t1);
+
+ if( !--k )
+ {
+ k = 3;
+ p -= 12;
+ }
+ }
+ }
+ else
+ {
+ for( i = 0; i <= size.width - 4; i += 4 )
+ {
+ int t0, t1, t;
+
+ temp = RNG_NEXT(temp);
+ t = (int)temp;
+ t0 = (t & p[i + 12]) + p[i];
+ t1 = ((t >> 8) & p[i + 13]) + p[i+1];
+ arr[i] = saturate_cast<T>(t0);
+ arr[i+1] = saturate_cast<T>(t1);
+
+ t0 = ((t >> 16) & p[i + 14]) + p[i + 2];
+ t1 = ((t >> 24) & p[i + 15]) + p[i + 3];
+ arr[i+2] = saturate_cast<T>(t0);
+ arr[i+3] = saturate_cast<T>(t1);
+
+ if( !--k )
+ {
+ k = 3;
+ p -= 12;
+ }
+ }
+ }
+
+ for( ; i < size.width; i++ )
+ {
+ unsigned t0;
+ temp = RNG_NEXT(temp);
+
+ t0 = ((int)temp & p[i + 12]) + p[i];
+ arr[i] = saturate_cast<T>(t0);
+ }
+ }
+
+ *state = temp;
+}
+
+
+template<typename T, typename PT> static void
+Randi_( Mat& _arr, uint64* state, const void* _param )
+{
+ uint64 temp = *state;
+ const PT* param = (const PT*)_param;
+ Size size = getContinuousSize(_arr,_arr.channels());
+
+ for( int y = 0; y < size.height; y++ )
+ {
+ T* arr = (T*)(_arr.data + _arr.step*y);
+ int i, k = 3;
+ const PT* p = param;
+
+ for( i = 0; i <= size.width - 4; i += 4 )
+ {
+ PT f0, f1;
+ temp = RNG_NEXT(temp);
+ f0 = (int)temp * p[i+12] + p[i];
+ temp = RNG_NEXT(temp);
+ f1 = (int)temp * p[i+13] + p[i+1];
+ arr[i] = saturate_cast<T>(cvFloor(f0));
+ arr[i+1] = saturate_cast<T>(cvFloor(f1));
+
+ temp = RNG_NEXT(temp);
+ f0 = (int)temp * p[i+14] + p[i+2];
+ temp = RNG_NEXT(temp);
+ f1 = (int)temp * p[i+15] + p[i+3];
+ arr[i+2] = saturate_cast<T>(cvFloor(f0));
+ arr[i+3] = saturate_cast<T>(cvFloor(f1));
+
+ if( !--k )
+ {
+ k = 3;
+ p -= 12;
+ }
+ }
+
+ for( ; i < size.width; i++ )
+ {
+ temp = RNG_NEXT(temp);
+ arr[i] = saturate_cast<T>(cvFloor((int)temp * p[i + 12] + p[i]));
+ }
+ }
+
+ *state = temp;
+}
+
+
+static void Randf_( Mat& _arr, uint64* state, const void* _param )
+{
+ uint64 temp = *state;
+ const float* param = (const float*)_param;
+ Size size = getContinuousSize(_arr,_arr.channels());
+
+ for( int y = 0; y < size.height; y++ )
+ {
+ float* arr = (float*)(_arr.data + _arr.step*y);
+ int i, k = 3;
+ const float* p = param;
+ for( i = 0; i <= size.width - 4; i += 4 )
+ {
+ float f0, f1;
+
+ temp = RNG_NEXT(temp);
+ f0 = (int)temp*p[i+12] + p[i];
+ temp = RNG_NEXT(temp);
+ f1 = (int)temp*p[i+13] + p[i+1];
+ arr[i] = f0; arr[i+1] = f1;
+
+ temp = RNG_NEXT(temp);
+ f0 = (int)temp*p[i+14] + p[i+2];
+ temp = RNG_NEXT(temp);
+ f1 = (int)temp*p[i+15] + p[i+3];
+ arr[i+2] = f0; arr[i+3] = f1;
+
+ if( !--k )
+ {
+ k = 3;
+ p -= 12;
+ }
+ }
+
+ for( ; i < size.width; i++ )
+ {
+ temp = RNG_NEXT(temp);
+ arr[i] = (int)temp*p[i+12] + p[i];
+ }
+ }
+
+ *state = temp;
+}
+
+
+static void
+Randd_( Mat& _arr, uint64* state, const void* _param )
+{
+ uint64 temp = *state;
+ const double* param = (const double*)_param;
+ Size size = getContinuousSize(_arr,_arr.channels());
+ int64 v = 0;
+
+ for( int y = 0; y < size.height; y++ )
+ {
+ double* arr = (double*)(_arr.data + _arr.step*y);
+ int i, k = 3;
+ const double* p = param;
+
+ for( i = 0; i <= size.width - 4; i += 4 )
+ {
+ double f0, f1;
+
+ temp = RNG_NEXT(temp);
+ v = (temp >> 32)|(temp << 32);
+ f0 = v*p[i+12] + p[i];
+ temp = RNG_NEXT(temp);
+ v = (temp >> 32)|(temp << 32);
+ f1 = v*p[i+12] + p[i];
+ arr[i] = f0; arr[i+1] = f1;
+
+ temp = RNG_NEXT(temp);
+ v = (temp >> 32)|(temp << 32);
+ f0 = v*p[i+12] + p[i];
+ temp = RNG_NEXT(temp);
+ v = (temp >> 32)|(temp << 32);
+ f1 = v*p[i+12] + p[i];
+ arr[i+2] = f0; arr[i+3] = f1;
+
+ if( !--k )
+ {
+ k = 3;
+ p -= 12;
+ }
+ }
+
+ for( ; i < size.width; i++ )
+ {
+ temp = RNG_NEXT(temp);
+ v = (temp >> 32)|(temp << 32);
+ arr[i] = v*p[i+12] + p[i];
+ }
+ }
+
+ *state = temp;
+}
+
+
+/*
+ The code below implements the algorithm described in
+ "The Ziggurat Method for Generating Random Variables"
+ by Marsaglia and Tsang, Journal of Statistical Software.
+*/
+static void
+Randn_0_1_32f_C1R( float* arr, int len, uint64* state )
+{
+ const float r = 3.442620f; // The start of the right tail
+ const float rng_flt = 2.3283064365386962890625e-10f; // 2^-32
+ static unsigned kn[127];
+ static float wn[128], fn[128];
+ uint64 temp = *state;
+ static bool initialized=false;
+ int i;
+
+ if( !initialized )
+ {
+ const double m1 = 2147483648.0;
+ double dn = 3.442619855899, tn = dn, vn = 9.91256303526217e-3;
+
+ // Set up the tables
+ double q = vn/std::exp(-.5*dn*dn);
+ kn[0] = (unsigned)((dn/q)*m1);
+ kn[1] = 0;
+
+ wn[0] = (float)(q/m1);
+ wn[127] = (float)(dn/m1);
+
+ fn[0] = 1.f;
+ fn[127] = (float)std::exp(-.5*dn*dn);
+
+ for(i=126;i>=1;i--)
+ {
+ dn = std::sqrt(-2.*std::log(vn/dn+std::exp(-.5*dn*dn)));
+ kn[i+1] = (unsigned)((dn/tn)*m1);
+ tn = dn;
+ fn[i] = (float)std::exp(-.5*dn*dn);
+ wn[i] = (float)(dn/m1);
+ }
+ initialized = true;
+ }
+
+ for( i = 0; i < len; i++ )
+ {
+ float x, y;
+ for(;;)
+ {
+ int hz = (int)temp;
+ temp = RNG_NEXT(temp);
+ int iz = hz & 127;
+ x = hz*wn[iz];
+ if( (unsigned)std::abs(hz) < kn[iz] )
+ break;
+ if( iz == 0) // iz==0, handles the base strip
+ {
+ do
+ {
+ x = (unsigned)temp*rng_flt;
+ temp = RNG_NEXT(temp);
+ y = (unsigned)temp*rng_flt;
+ temp = RNG_NEXT(temp);
+ x = (float)(-std::log(x+FLT_MIN)*0.2904764);
+ y = (float)-std::log(y+FLT_MIN);
+ } // .2904764 is 1/r
+ while( y + y < x*x );
+ x = hz > 0 ? r + x : -r - x;
+ break;
+ }
+ // iz > 0, handle the wedges of other strips
+ y = (unsigned)temp*rng_flt;
+ temp = RNG_NEXT(temp);
+ if( fn[iz] + y*(fn[iz - 1] - fn[iz]) < std::exp(-.5*x*x) )
+ break;
+ }
+ arr[i] = x;
+ }
+ *state = temp;
+}
+
+
+template<typename T, typename PT> static void
+Randn_( Mat& _arr, uint64* state, const void* _param )
+{
+ const int RAND_BUF_SIZE = 96;
+ float buffer[RAND_BUF_SIZE];
+ const PT* param = (const PT*)_param;
+ Size size = getContinuousSize(_arr, _arr.channels());
+
+ for( int y = 0; y < size.height; y++ )
+ {
+ T* arr = (T*)(_arr.data + _arr.step*y);
+ int i, j, len = RAND_BUF_SIZE;
+ for( i = 0; i < size.width; i += RAND_BUF_SIZE )
+ {
+ int k = 3;
+ const PT* p = param;
+
+ if( i + len > size.width )
+ len = size.width - i;
+
+ Randn_0_1_32f_C1R( buffer, len, state );
+
+ for( j = 0; j <= len - 4; j += 4 )
+ {
+ PT f0, f1;
+
+ f0 = buffer[j]*p[j+12] + p[j];
+ f1 = buffer[j+1]*p[j+13] + p[j+1];
+ arr[i+j] = saturate_cast<T>(f0);
+ arr[i+j+1] = saturate_cast<T>(f1);
+
+ f0 = buffer[j+2]*p[j+14] + p[j+2];
+ f1 = buffer[j+3]*p[j+15] + p[j+3];
+ arr[i+j+2] = saturate_cast<T>(f0);
+ arr[i+j+3] = saturate_cast<T>(f1);
+
+ if( --k == 0 )
+ {
+ k = 3;
+ p -= 12;
+ }
+ }
+
+ for( ; j < len; j++ )
+ arr[i+j] = saturate_cast<T>(buffer[j]*p[j+12] + p[j]);
+ }
+ }
+}
+
+
+typedef void (*RandFunc)(Mat& dst, uint64* state, const void* param);
+
+void RNG::fill( Mat& mat, int disttype, const Scalar& param1, const Scalar& param2 )
+{
+ static RandFunc rngtab[3][8] =
+ {
+ {
+ RandBits_<uchar>,
+ RandBits_<schar>,
+ RandBits_<ushort>,
+ RandBits_<short>,
+ RandBits_<int>, 0, 0, 0},
+
+ {Randi_<uchar,float>,
+ Randi_<schar,float>,
+ Randi_<ushort,float>,
+ Randi_<short,float>,
+ Randi_<int,float>,
+ Randf_, Randd_, 0},
+
+ {Randn_<uchar,float>,
+ Randn_<schar,float>,
+ Randn_<ushort,float>,
+ Randn_<short,float>,
+ Randn_<int,float>,
+ Randn_<float,float>,
+ Randn_<double,double>, 0}
+ };
+
+ int depth = mat.depth(), channels = mat.channels();
+ double dparam[2][12];
+ float fparam[2][12];
+ int iparam[2][12];
+ void* param = dparam;
+ int i, fast_int_mode = 0;
+ RandFunc func = 0;
+
+ CV_Assert( channels <= 4 );
+
+ if( disttype == UNIFORM )
+ {
+ if( depth <= CV_32S )
+ {
+ for( i = 0, fast_int_mode = 1; i < channels; i++ )
+ {
+ int t0 = iparam[0][i] = cvCeil(param1.val[i]);
+ int t1 = iparam[1][i] = cvFloor(param2.val[i]) - t0;
+ double diff = param1.val[i] - param2.val[i];
+
+ fast_int_mode &= INT_MIN <= diff && diff <= INT_MAX && (t1 & (t1 - 1)) == 0;
+ }
+ }
+
+ if( fast_int_mode )
+ {
+ for( i = 0; i < channels; i++ )
+ iparam[1][i]--;
+
+ for( ; i < 12; i++ )
+ {
+ int t0 = iparam[0][i - channels];
+ int t1 = iparam[1][i - channels];
+
+ iparam[0][i] = t0;
+ iparam[1][i] = t1;
+ }
+
+ func = rngtab[0][depth];
+ param = iparam;
+ }
+ else
+ {
+ double scale = depth == CV_64F ?
+ 5.4210108624275221700372640043497e-20 : // 2**-64
+ 2.3283064365386962890625e-10; // 2**-32
+
+ // for each channel i compute such dparam[0][i] & dparam[1][i],
+ // so that a signed 32/64-bit integer X is transformed to
+ // the range [param1.val[i], param2.val[i]) using
+ // dparam[1][i]*X + dparam[0][i]
+ for( i = 0; i < channels; i++ )
+ {
+ double t0 = param1.val[i];
+ double t1 = param2.val[i];
+ dparam[0][i] = (t1 + t0)*0.5;
+ dparam[1][i] = (t1 - t0)*scale;
+ }
+
+ func = rngtab[1][depth];
+ param = dparam;
+ }
+ }
+ else if( disttype == CV_RAND_NORMAL )
+ {
+ for( i = 0; i < channels; i++ )
+ {
+ double t0 = param1.val[i];
+ double t1 = param2.val[i];
+
+ dparam[0][i] = t0;
+ dparam[1][i] = t1;
+ }
+
+ func = rngtab[2][depth];
+ param = dparam;
+ }
+ else
+ CV_Error( CV_StsBadArg, "Unknown distribution type" );
+
+ if( !fast_int_mode )
+ {
+ for( i = channels; i < 12; i++ )
+ {
+ double t0 = dparam[0][i - channels];
+ double t1 = dparam[1][i - channels];
+
+ dparam[0][i] = t0;
+ dparam[1][i] = t1;
+ }
+
+ if( depth != CV_64F )
+ {
+ for( i = 0; i < 12; i++ )
+ {
+ fparam[0][i] = (float)dparam[0][i];
+ fparam[1][i] = (float)dparam[1][i];
+ }
+ param = fparam;
+ }
+ }
+
+ CV_Assert( func != 0);
+
+ func( mat, &state, param );
+}
+
+
+#ifdef WIN32
+#ifdef WINCE
+# define TLS_OUT_OF_INDEXES ((DWORD)0xFFFFFFFF)
+#endif
+static DWORD tlsRNGKey = TLS_OUT_OF_INDEXES;
+
+void deleteThreadRNGData()
+{
+ if( tlsRNGKey != TLS_OUT_OF_INDEXES )
+ delete (RNG*)TlsGetValue( tlsRNGKey );
+}
+
+RNG& theRNG()
+{
+ if( tlsRNGKey == TLS_OUT_OF_INDEXES )
+ {
+ tlsRNGKey = TlsAlloc();
+ CV_Assert(tlsRNGKey != TLS_OUT_OF_INDEXES);
+ }
+ RNG* rng = (RNG*)TlsGetValue( tlsRNGKey );
+ if( !rng )
+ {
+ rng = new RNG;
+ TlsSetValue( tlsRNGKey, rng );
+ }
+ return *rng;
+}
+
+#else
+
+static pthread_key_t tlsRNGKey = 0;
+
+static void deleteRNG(void* data)
+{
+ delete (RNG*)data;
+}
+
+RNG& theRNG()
+{
+ if( !tlsRNGKey )
+ {
+ pthread_key_create(&tlsRNGKey, deleteRNG);
+ CV_Assert(tlsRNGKey != 0);
+ }
+ RNG* rng = (RNG*)pthread_getspecific(tlsRNGKey);
+ if( !rng )
+ {
+ rng = new RNG;
+ pthread_setspecific(tlsRNGKey, rng);
+ }
+ return *rng;
+}
+
+#endif
+
+template<typename T> static void
+randShuffle_( Mat& _arr, RNG& rng, double iterFactor )
+{
+ int sz = _arr.rows*_arr.cols, iters = cvRound(iterFactor*sz);
+ if( _arr.isContinuous() )
+ {
+ T* arr = (T*)_arr.data;
+ for( int i = 0; i < iters; i++ )
+ {
+ int j = (unsigned)rng % sz, k = (unsigned)rng % sz;
+ std::swap( arr[j], arr[k] );
+ }
+ }
+ else
+ {
+ uchar* data = _arr.data;
+ size_t step = _arr.step;
+ int cols = _arr.cols;
+ for( int i = 0; i < iters; i++ )
+ {
+ int j1 = (unsigned)rng % sz, k1 = (unsigned)rng % sz;
+ int j0 = j1/cols, k0 = k1/cols;
+ j1 -= j0*cols; k1 -= k0*cols;
+ std::swap( ((T*)(data + step*j0))[j1], ((T*)(data + step*k0))[k1] );
+ }
+ }
+}
+
+typedef void (*RandShuffleFunc)( Mat& dst, RNG& rng, double iterFactor );
+
+void randShuffle( Mat& dst, double iterFactor, RNG* _rng )
+{
+ RandShuffleFunc tab[] =
+ {
+ 0,
+ randShuffle_<uchar>, // 1
+ randShuffle_<ushort>, // 2
+ randShuffle_<Vec<uchar,3> >, // 3
+ randShuffle_<int>, // 4
+ 0,
+ randShuffle_<Vec<ushort,3> >, // 6
+ 0,
+ randShuffle_<int64>, // 8
+ 0, 0, 0,
+ randShuffle_<Vec<int,3> >, // 12
+ 0, 0, 0,
+ randShuffle_<Vec<int64,2> >, // 16
+ 0, 0, 0, 0, 0, 0, 0,
+ randShuffle_<Vec<int64,3> >, // 24
+ 0, 0, 0, 0, 0, 0, 0,
+ randShuffle_<Vec<int64,4> > // 32
+ };
+
+ RNG& rng = _rng ? *_rng : theRNG();
+ CV_Assert( dst.elemSize() <= 32 );
+ RandShuffleFunc func = tab[dst.elemSize()];
+ CV_Assert( func != 0 );
+ func( dst, rng, iterFactor );
+}
+
+}
+
+CV_IMPL void
+cvRandArr( CvRNG* _rng, CvArr* arr, int disttype, CvScalar param1, CvScalar param2 )
+{
+ cv::Mat mat = cv::cvarrToMat(arr);
+ // !!! this will only work for current 64-bit MWC RNG !!!
+ cv::RNG& rng = _rng ? (cv::RNG&)*_rng : cv::theRNG();
+ rng.fill(mat, disttype == CV_RAND_NORMAL ?
+ cv::RNG::NORMAL : cv::RNG::UNIFORM, param1, param2 );
+}
+
+CV_IMPL void cvRandShuffle( CvArr* arr, CvRNG* _rng, double iter_factor )
+{
+ cv::Mat dst = cv::cvarrToMat(arr);
+ cv::RNG& rng = _rng ? (cv::RNG&)*_rng : cv::theRNG();
+ cv::randShuffle( dst, iter_factor, &rng );
+}
+
+/* End of file. */
projects / opencv / blobdiff