From f454929d9c2f96f74cee0c8006458fd40e20da46 Mon Sep 17 00:00:00 2001 From: Olexa Bilaniuk Date: Sat, 21 Feb 2015 12:31:55 -0500 Subject: [PATCH] PRNG changes: xorshift128+ algorithm, and seeding API. - Switched to the extremely fast, while simple and high-quality, xorshift128+ PRNG algorithm by Sebastiano Vigna in "Further scramblings of Marsaglia's xorshift generators. CoRR, abs/1402.6246, 2014" (2^128-1 period, passes BigCrush tests). Performance improved by 10% over random(). - Added an API to allow seeding with a specified seed, rather than using rand() or random(). This allows deterministic, reproducible results in tests using our algorithm (although findHomography() does not yet support passing an entropy source on its own end). --- modules/calib3d/src/rhorefc.cpp | 214 ++++++++++++++++++++------------ modules/calib3d/src/rhorefc.h | 14 +++ 2 files changed, 149 insertions(+), 79 deletions(-) diff --git a/modules/calib3d/src/rhorefc.cpp b/modules/calib3d/src/rhorefc.cpp index f60f68a04..86d2f0977 100644 --- a/modules/calib3d/src/rhorefc.cpp +++ b/modules/calib3d/src/rhorefc.cpp @@ -161,6 +161,11 @@ struct RHO_HEST_REFC{ float* Jte; /* Jte vector */ } lm; + /* PRNG XORshift128+ */ + struct{ + uint64_t s[2]; /* PRNG state */ + } prng; + /* Initialized? */ int init; @@ -205,6 +210,11 @@ struct RHO_HEST_REFC{ inline int PROSACPhaseEndReached(void); inline void PROSACGoToNextPhase(void); inline void getPROSACSample(void); + inline void rndSmpl(unsigned sampleSize, + unsigned* currentSample, + unsigned dataSetSize); + inline double fastRandom(void); + inline void fastSeed(uint64_t seed); inline int isSampleDegenerate(void); inline void generateModel(void); inline int isModelDegenerate(void); @@ -239,10 +249,6 @@ static inline void sacInitNonRand (double beta, static inline double sacInitPEndFpI (const unsigned ransacConvg, const unsigned n, const unsigned s); -static inline void sacRndSmpl (unsigned sampleSize, - unsigned* currentSample, - unsigned dataSetSize); -static inline double sacRandom (void); static inline unsigned sacCalcIterBound (double confidence, double inlierRate, unsigned sampleSize, @@ -304,6 +310,15 @@ int rhoRefCEnsureCapacity(RHO_HEST_REFC* p, unsigned N, double beta){ } +/** + * Seeds the internal PRNG with the given seed. + */ + +void rhoRefCSeed(RHO_HEST_REFC* p, unsigned long long seed){ + p->fastSeed((uint64_t)seed); +} + + /** * External access to context destructor. * @@ -459,6 +474,12 @@ inline int RHO_HEST_REFC::initialize(void){ lm.tmp1 = NULL; lm.Jte = NULL; +#ifdef _WIN32 + fastSeed(rand()); +#else + fastSeed(random()); +#endif + int areAllAllocsSuccessful = ctrl.smpl && curr.H && @@ -950,13 +971,121 @@ inline void RHO_HEST_REFC::PROSACGoToNextPhase(void){ inline void RHO_HEST_REFC::getPROSACSample(void){ if(ctrl.i > ctrl.phEndI){ /* FIXME: Dubious. Review. */ - sacRndSmpl(4, ctrl.smpl, ctrl.phNum);/* Used to be phMax */ + rndSmpl(4, ctrl.smpl, ctrl.phNum);/* Used to be phMax */ }else{ - sacRndSmpl(3, ctrl.smpl, ctrl.phNum-1); + rndSmpl(3, ctrl.smpl, ctrl.phNum-1); ctrl.smpl[3] = ctrl.phNum-1; } } +/** + * Choose, without repetition, sampleSize integers in the range [0, numDataPoints). + */ + +inline void RHO_HEST_REFC::rndSmpl(unsigned sampleSize, + unsigned* currentSample, + unsigned dataSetSize){ + /** + * If sampleSize is very close to dataSetSize, we use selection sampling. + * Otherwise we use the naive sampling technique wherein we select random + * indexes until sampleSize of them are distinct. + */ + + if(sampleSize*2>dataSetSize){ + /** + * Selection Sampling: + * + * Algorithm S (Selection sampling technique). To select n records at random from a set of N, where 0 < n ≤ N. + * S1. [Initialize.] Set t ← 0, m ← 0. (During this algorithm, m represents the number of records selected so far, + * and t is the total number of input records that we have dealt with.) + * S2. [Generate U.] Generate a random number U, uniformly distributed between zero and one. + * S3. [Test.] If (N – t)U ≥ n – m, go to step S5. + * S4. [Select.] Select the next record for the sample, and increase m and t by 1. If m < n, go to step S2; + * otherwise the sample is complete and the algorithm terminates. + * S5. [Skip.] Skip the next record (do not include it in the sample), increase t by 1, and go back to step S2. + * + * Replaced m with i and t with j in the below code. + */ + + unsigned i=0,j=0; + + for(i=0;i> 17; // b + x ^= y ^ (y >> 26); // c + prng.s[0] = y; + prng.s[1] = x; + uint64_t s = x + y; + + return s * 5.421010862427522e-20;/* 2^-64 */ +} + +/** + * Seeds the PRNG. + * + * The seed should not be zero, since the state must be initialized to non-zero. + */ + +inline void RHO_HEST_REFC::fastSeed(uint64_t seed){ + int i; + + prng.s[0] = seed; + prng.s[1] = ~seed;/* Guarantees one of the elements will be non-zero. */ + + /** + * Escape from zero-land (see xorshift128+ paper). Approximately 20 + * iterations required according to the graph. + */ + + for(i=0;i<20;i++){ + fastRandom(); + } +} + /** * Checks whether the *sample* is degenerate prior to model generation. * - First, the extremely cheap numerical degeneracy test is run, which weeds @@ -1395,79 +1524,6 @@ static inline double sacInitPEndFpI(const unsigned ransacConvg, return ransacConvg*numer/denom; } -/** - * Choose, without repetition, sampleSize integers in the range [0, numDataPoints). - */ - -static inline void sacRndSmpl(unsigned sampleSize, - unsigned* currentSample, - unsigned dataSetSize){ - /** - * If sampleSize is very close to dataSetSize, we use selection sampling. - * Otherwise we use the naive sampling technique wherein we select random - * indexes until sampleSize of them are distinct. - */ - - if(sampleSize*2>dataSetSize){ - /** - * Selection Sampling: - * - * Algorithm S (Selection sampling technique). To select n records at random from a set of N, where 0 < n ≤ N. - * S1. [Initialize.] Set t ← 0, m ← 0. (During this algorithm, m represents the number of records selected so far, - * and t is the total number of input records that we have dealt with.) - * S2. [Generate U.] Generate a random number U, uniformly distributed between zero and one. - * S3. [Test.] If (N – t)U ≥ n – m, go to step S5. - * S4. [Select.] Select the next record for the sample, and increase m and t by 1. If m < n, go to step S2; - * otherwise the sample is complete and the algorithm terminates. - * S5. [Skip.] Skip the next record (do not include it in the sample), increase t by 1, and go back to step S2. - * - * Replaced m with i and t with j in the below code. - */ - - unsigned i=0,j=0; - - for(i=0;i