[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

# vigra/random.hxx

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*                  Copyright 2008 by Ullrich Koethe                    */
    00004 /*                                                                      */
    00005 /*    This file is part of the VIGRA computer vision library.           */
    00006 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
    00007 /*    The VIGRA Website is                                              */
    00008 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00009 /*    Please direct questions, bug reports, and contributions to        */
    00010 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00011 /*        vigra@informatik.uni-hamburg.de                               */
    00012 /*                                                                      */
    00013 /*    Permission is hereby granted, free of charge, to any person       */
    00014 /*    obtaining a copy of this software and associated documentation    */
    00015 /*    files (the "Software"), to deal in the Software without           */
    00016 /*    restriction, including without limitation the rights to use,      */
    00017 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00018 /*    sell copies of the Software, and to permit persons to whom the    */
    00019 /*    Software is furnished to do so, subject to the following          */
    00020 /*    conditions:                                                       */
    00021 /*                                                                      */
    00022 /*    The above copyright notice and this permission notice shall be    */
    00023 /*    included in all copies or substantial portions of the             */
    00024 /*    Software.                                                         */
    00025 /*                                                                      */
    00026 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00027 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00028 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00029 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00030 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00031 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00032 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00033 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00034 /*                                                                      */
    00035 /************************************************************************/
    00036
    00037
    00038 #ifndef VIGRA_RANDOM_HXX
    00039 #define VIGRA_RANDOM_HXX
    00040
    00041 #include "mathutil.hxx"
    00042 #include "functortraits.hxx"
    00043
    00044 #include <time.h>
    00045
    00046 namespace vigra {
    00047
    00048 enum RandomSeedTag { RandomSeed };
    00049
    00050 namespace detail {
    00051
    00052 enum RandomEngineTag { TT800, MT19937 };
    00053
    00054 template<RandomEngineTag EngineTag>
    00055 struct RandomState;
    00056
    00057 template <RandomEngineTag EngineTag>
    00058 void seed(UInt32 theSeed, RandomState<EngineTag> & engine)
    00059 {
    00060     engine.state_[0] = theSeed;
    00061     for(UInt32 i=1; i<RandomState<EngineTag>::N; ++i)
    00062     {
    00063         engine.state_[i] = 1812433253U * (engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) + i;
    00064     }
    00065 }
    00066
    00067 template <class Iterator, RandomEngineTag EngineTag>
    00068 void seed(Iterator init, UInt32 key_length, RandomState<EngineTag> & engine)
    00069 {
    00070     const UInt32 N = RandomState<EngineTag>::N;
    00071     int k = (int)std::max(N, key_length);
    00072     UInt32 i = 1, j = 0;
    00073     Iterator data = init;
    00074     for (; k; --k)
    00075     {
    00076         engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1664525U))
    00077                            + *data + j; /* non linear */
    00078         ++i; ++j; ++data;
    00079
    00080         if (i >= N)
    00081         {
    00082             engine.state_[0] = engine.state_[N-1];
    00083             i=1;
    00084         }
    00085         if (j>=key_length)
    00086         {
    00087             j=0;
    00088             data = init;
    00089         }
    00090     }
    00091
    00092     for (k=N-1; k; --k)
    00093     {
    00094         engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1566083941U))
    00095                            - i; /* non linear */
    00096         ++i;
    00097         if (i>=N)
    00098         {
    00099             engine.state_[0] = engine.state_[N-1];
    00100             i=1;
    00101         }
    00102     }
    00103
    00104     engine.state_[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
    00105 }
    00106
    00107 template <RandomEngineTag EngineTag>
    00108 void seed(RandomSeedTag, RandomState<EngineTag> & engine)
    00109 {
    00110     UInt32 init[2] = { (UInt32)time(0), (UInt32)clock() };
    00111     seed(init, 2, engine);
    00112 }
    00113
    00114
    00115     /* Tempered twister TT800 by M. Matsumoto */
    00116 template<>
    00117 struct RandomState<TT800>
    00118 {
    00119     static const UInt32 N = 25, M = 7;
    00120
    00121     mutable UInt32 state_[N];
    00122     mutable UInt32 current_;
    00123
    00124     RandomState()
    00125     : current_(0)
    00126     {
    00127         UInt32 seeds[N] = {
    00128             0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23,
    00129             0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825,
    00130             0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f,
    00131             0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9,
    00132             0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb
    00133         };
    00134
    00135         for(UInt32 i=0; i<N; ++i)
    00136             state_[i] = seeds[i];
    00137     }
    00138
    00139   protected:
    00140
    00141     UInt32 get() const
    00142     {
    00143         if(current_ == N)
    00144             generateNumbers();
    00145
    00146         UInt32 y = state_[current_++];
    00147         y ^= (y << 7) & 0x2b5b2500;
    00148         y ^= (y << 15) & 0xdb8b0000;
    00149         return y ^ (y >> 16);
    00150     }
    00151
    00152     void generateNumbers() const;
    00153
    00154     void seedImpl(RandomSeedTag)
    00155     {
    00156         seed(RandomSeed, *this);
    00157     }
    00158
    00159     void seedImpl(UInt32 theSeed)
    00160     {
    00161         seed(theSeed, *this);
    00162     }
    00163
    00164     template<class Iterator>
    00165     void seedImpl(Iterator init, UInt32 length)
    00166     {
    00167         seed(init, length, *this);
    00168     }
    00169 };
    00170
    00171 void RandomState<TT800>::generateNumbers() const
    00172 {
    00173     UInt32 mag01[2]= { 0x0, 0x8ebfd028 };
    00174
    00175     for(UInt32 i=0; i<N-M; ++i)
    00176     {
    00177         state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
    00178     }
    00179     for (UInt32 i=N-M; i<N; ++i)
    00180     {
    00181         state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
    00182     }
    00183     current_ = 0;
    00184 }
    00185
    00186     /* Mersenne twister MT19937 by M. Matsumoto */
    00187 template<>
    00188 struct RandomState<MT19937>
    00189 {
    00190     static const UInt32 N = 624, M = 397;
    00191
    00192     mutable UInt32 state_[N];
    00193     mutable UInt32 current_;
    00194
    00195     RandomState()
    00196     : current_(0)
    00197     {
    00198         seed(19650218U, *this);
    00199     }
    00200
    00201   protected:
    00202
    00203     UInt32 get() const
    00204     {
    00205         if(current_ == N)
    00206             generateNumbers();
    00207
    00208         UInt32 x = state_[current_++];
    00209         x ^= (x >> 11);
    00210         x ^= (x << 7) & 0x9D2C5680U;
    00211         x ^= (x << 15) & 0xEFC60000U;
    00212         return x ^ (x >> 18);
    00213     }
    00214
    00215     void generateNumbers() const;
    00216
    00217     static UInt32 twiddle(UInt32 u, UInt32 v)
    00218     {
    00219         return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
    00220                 ^ ((v & 1U) ? 0x9908B0DFU : 0x0U);
    00221     }
    00222
    00223     void seedImpl(RandomSeedTag)
    00224     {
    00225         seed(RandomSeed, *this);
    00226         generateNumbers();
    00227     }
    00228
    00229     void seedImpl(UInt32 theSeed)
    00230     {
    00231         seed(theSeed, *this);
    00232         generateNumbers();
    00233     }
    00234
    00235     template<class Iterator>
    00236     void seedImpl(Iterator init, UInt32 length)
    00237     {
    00238         seed(19650218U, *this);
    00239         seed(init, length, *this);
    00240         generateNumbers();
    00241     }
    00242 };
    00243
    00244 void RandomState<MT19937>::generateNumbers() const
    00245 {
    00246     for (int i = 0; i < (N - M); ++i)
    00247     {
    00248         state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]);
    00249     }
    00250     for (int i = N - M; i < (N - 1); ++i)
    00251     {
    00252         state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]);
    00253     }
    00254     state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]);
    00255     current_ = 0;
    00256 }
    00257
    00258 } // namespace detail
    00259
    00260
    00261 /** \addtogroup RandomNumberGeneration Random Number Generation
    00262
    00263      High-quality random number generators and related functors.
    00264 */
    00265 //@{
    00266
    00267 /** Generic random number generator.
    00268
    00269     The actual generator is passed in the template argument <tt>Engine</tt>. Two generators
    00270     are currently available:
    00271     <ul>
    00272     <li> <tt>RandomMT19937</tt>: The state-of-the-art <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">Mersenne Twister</a> with a state length of 2<sup>19937</sup> and very high statistical quality.
    00273     <li> <tt>RandomTT800</tt>: (default) The Tempered Twister, a simpler predecessor of the Mersenne Twister with period length 2<sup>800</sup>.
    00274     </ul>
    00275
    00276     Both generators have been designed by <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/eindex.html">Makoto Matsumoto</a>.
    00277
    00278     <b>Traits defined:</b>
    00279
    00280     <tt>FunctorTraits<RandomNumberGenerator<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
    00281 */
    00282 template <class Engine = detail::RandomState<detail::TT800> >
    00283 class RandomNumberGenerator
    00284 : public Engine
    00285 {
    00286     mutable double normalCached_;
    00287     mutable bool normalCachedValid_;
    00288
    00289   public:
    00290
    00291         /** Create a new random generator object with standard seed.
    00292
    00293             Due to standard seeding, the random numbers generated will always be the same.
    00294             This is useful for debugging.
    00295         */
    00296     RandomNumberGenerator()
    00297     : normalCached_(0.0),
    00298       normalCachedValid_(false)
    00299     {}
    00300
    00301         /** Create a new random generator object with a random seed.
    00302
    00303             The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
    00304             values.
    00305
    00306             <b>Usage:</b>
    00307             \code
    00308             RandomNumberGenerator<> rnd = RandomNumberGenerator<>(RandomSeed);
    00309             \endcode
    00310         */
    00311     RandomNumberGenerator(RandomSeedTag)
    00312     : normalCached_(0.0),
    00313       normalCachedValid_(false)
    00314     {
    00315         this->seedImpl(RandomSeed);
    00316     }
    00317
    00318         /** Create a new random generator object from the given seed.
    00319
    00320             The same seed will always produce identical random sequences.
    00321         */
    00322     RandomNumberGenerator(UInt32 theSeed)
    00323     : normalCached_(0.0),
    00324       normalCachedValid_(false)
    00325     {
    00326         this->seedImpl(theSeed);
    00327     }
    00328
    00329         /** Create a new random generator object from the given seed sequence.
    00330
    00331             Longer seed sequences lead to better initialization in the sense that the generator's
    00332             state space is covered much better than is possible with 32-bit seeds alone.
    00333         */
    00334     template<class Iterator>
    00335     RandomNumberGenerator(Iterator init, UInt32 length)
    00336     : normalCached_(0.0),
    00337       normalCachedValid_(false)
    00338     {
    00339         this->seedImpl(init, length);
    00340     }
    00341
    00342
    00343         /** Re-initialize the random generator object with a random seed.
    00344
    00345             The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
    00346             values.
    00347
    00348             <b>Usage:</b>
    00349             \code
    00350             RandomNumberGenerator<> rnd = ...;
    00351             ...
    00352             rnd.seed(RandomSeed);
    00353             \endcode
    00354         */
    00355     void seed(RandomSeedTag)
    00356     {
    00357         this->seedImpl(RandomSeed);
    00358         normalCachedValid_ = false;
    00359     }
    00360
    00361         /** Re-initialize the random generator object from the given seed.
    00362
    00363             The same seed will always produce identical random sequences.
    00364         */
    00365     void seed(UInt32 theSeed)
    00366     {
    00367         this->seedImpl(theSeed);
    00368         normalCachedValid_ = false;
    00369     }
    00370
    00371         /** Re-initialize the random generator object from the given seed sequence.
    00372
    00373             Longer seed sequences lead to better initialization in the sense that the generator's
    00374             state space is covered much better than is possible with 32-bit seeds alone.
    00375         */
    00376     template<class Iterator>
    00377     void seed(Iterator init, UInt32 length)
    00378     {
    00379         this->seedImpl(init, length);
    00380         normalCachedValid_ = false;
    00381     }
    00382
    00383         /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
    00384
    00385             That is, 0 &lt;= i &lt; 2<sup>32</sup>.
    00386         */
    00387     UInt32 operator()() const
    00388     {
    00389         return this->get();
    00390     }
    00391
    00392         /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
    00393
    00394             That is, 0 &lt;= i &lt; 2<sup>32</sup>.
    00395         */
    00396     UInt32 uniformInt() const
    00397     {
    00398         return this->get();
    00399     }
    00400
    00401
    00402 #if 0 // difficult implementation necessary if low bits are not sufficiently random
    00403         // in [0,beyond)
    00404     UInt32 uniformInt(UInt32 beyond) const
    00405     {
    00406         if(beyond < 2)
    00407             return 0;
    00408
    00409         UInt32 factor = factorForUniformInt(beyond);
    00410         UInt32 res = this->get() / factor;
    00411
    00412         // Use rejection method to avoid quantization bias.
    00413         // On average, we will need two raw random numbers to generate one.
    00414         while(res >= beyond)
    00415             res = this->get() / factor;
    00416         return res;
    00417     }
    00418 #endif /* #if 0 */
    00419
    00420         /** Return a uniformly distributed integer random number in [0, <tt>beyond</tt>).
    00421
    00422             That is, 0 &lt;= i &lt; <tt>beyond</tt>.
    00423         */
    00424     UInt32 uniformInt(UInt32 beyond) const
    00425     {
    00426         // in [0,beyond) -- simple implementation when low bits are sufficiently random, which is
    00427         // the case for TT800 and MT19937
    00428         if(beyond < 2)
    00429             return 0;
    00430
    00431         UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond;
    00432         UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder;
    00433         UInt32 res = this->get();
    00434
    00435         // Use rejection method to avoid quantization bias.
    00436         // We will need two raw random numbers in amortized worst case.
    00437         while(res > lastSafeValue)
    00438             res = this->get();
    00439         return res % beyond;
    00440     }
    00441
    00442         /** Return a uniformly distributed double-precision random number in [0.0, 1.0).
    00443
    00444             That is, 0.0 &lt;= i &lt; 1.0. All 53-bit bits of the mantissa are random (two 32-bit integers are used to
    00445             create this number).
    00446         */
    00447     double uniform53() const
    00448     {
    00449         // make full use of the entire 53-bit mantissa of a double, by Isaku Wada
    00450         return ( (this->get() >> 5) * 67108864.0 + (this->get() >> 6)) * (1.0/9007199254740992.0);
    00451     }
    00452
    00453         /** Return a uniformly distributed double-precision random number in [0.0, 1.0].
    00454
    00455             That is, 0.0 &lt;= i &lt;= 1.0. This nuber is computed by <tt>uniformInt()</tt> / 2<sup>32</sup>,
    00456             so it has effectively only 32 random bits.
    00457         */
    00458     double uniform() const
    00459     {
    00460         return (double)this->get() / 4294967295.0;
    00461     }
    00462
    00463         /** Return a uniformly distributed double-precision random number in [lower, upper].
    00464
    00465             That is, <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>. This number is computed
    00466             from <tt>uniform()</tt>, so it has effectively only 32 random bits.
    00467         */
    00468     double uniform(double lower, double upper) const
    00469     {
    00470         vigra_precondition(lower < upper,
    00471           "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound.");
    00472         return uniform() * (upper-lower) + lower;
    00473     }
    00474
    00475         /** Return a standard normal variate (Gaussian) random number.
    00476
    00477             Mean is zero, standard deviation is 1.0. It uses the polar form of the
    00478             Box-Muller transform.
    00479         */
    00480     double normal() const;
    00481
    00482         /** Return a normal variate (Gaussian) random number with the given mean and standard deviation.
    00483
    00484             It uses the polar form of the Box-Muller transform.
    00485         */
    00486     double normal(double mean, double stddev) const
    00487     {
    00488         vigra_precondition(stddev > 0.0,
    00489           "RandomNumberGenerator::normal(): standard deviation must be positive.");
    00490         return normal()*stddev + mean;
    00491     }
    00492
    00493         /** Access the global (program-wide) instance of the present random number generator.
    00494
    00495             Normally, you will create a local generator by one of the constructor calls. But sometimes
    00496             it is useful to have all program parts access the same generator.
    00497         */
    00498     static RandomNumberGenerator & global()
    00499     {
    00500         static RandomNumberGenerator generator;
    00501         return generator;
    00502     }
    00503
    00504     static UInt32 factorForUniformInt(UInt32 range)
    00505     {
    00506         return (range > 2147483648U || range == 0)
    00507                      ? 1
    00508                      : 2*(2147483648U / ceilPower2(range));
    00509     }
    00510 };
    00511
    00512 template <class Engine>
    00513 double RandomNumberGenerator<Engine>::normal() const
    00514 {
    00515     if(normalCachedValid_)
    00516     {
    00517         normalCachedValid_ = false;
    00518         return normalCached_;
    00519     }
    00520     else
    00521     {
    00522         double x1, x2, w;
    00523         do
    00524         {
    00525              x1 = uniform(-1.0, 1.0);
    00526              x2 = uniform(-1.0, 1.0);
    00527              w = x1 * x1 + x2 * x2;
    00528         }
    00529         while ( w > 1.0 || w == 0.0);
    00530
    00531         w = std::sqrt( -2.0 * std::log( w )  / w );
    00532
    00533         normalCached_ = x2 * w;
    00534         normalCachedValid_ = true;
    00535
    00536         return x1 * w;
    00537     }
    00538 }
    00539
    00540     /** Shorthand for the TT800 random number generator class.
    00541     */
    00542 typedef RandomNumberGenerator<>  RandomTT800;
    00543
    00544     /** Shorthand for the MT19937 random number generator class.
    00545     */
    00546 typedef RandomNumberGenerator<detail::RandomState<detail::MT19937> > RandomMT19937;
    00547
    00548     /** Access the global (program-wide) instance of the TT800 random number generator.
    00549     */
    00550 inline RandomTT800   & randomTT800()   { return RandomTT800::global(); }
    00551
    00552     /** Access the global (program-wide) instance of the MT19937 random number generator.
    00553     */
    00554 inline RandomMT19937 & randomMT19937() { return RandomMT19937::global(); }
    00555
    00556 template <class Engine>
    00557 class FunctorTraits<RandomNumberGenerator<Engine> >
    00558 {
    00559   public:
    00560     typedef RandomNumberGenerator<Engine> type;
    00561
    00562     typedef VigraTrueType  isInitializer;
    00563
    00564     typedef VigraFalseType isUnaryFunctor;
    00565     typedef VigraFalseType isBinaryFunctor;
    00566     typedef VigraFalseType isTernaryFunctor;
    00567
    00568     typedef VigraFalseType isUnaryAnalyser;
    00569     typedef VigraFalseType isBinaryAnalyser;
    00570     typedef VigraFalseType isTernaryAnalyser;
    00571 };
    00572
    00573
    00574 /** Functor to create uniformely distributed integer random numbers.
    00575
    00576     This functor encapsulates the appropriate functions of the given random number
    00577     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
    00578     in an STL-compatible interface.
    00579
    00580     <b>Traits defined:</b>
    00581
    00582     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> and
    00583     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isUnaryFunctor</tt> are true (<tt>VigraTrueType</tt>).
    00584 */
    00585 template <class Engine = RandomTT800>
    00586 class UniformIntRandomFunctor
    00587 {
    00588     UInt32 lower_, difference_, factor_;
    00589     Engine & generator_;
    00590     bool useLowBits_;
    00591
    00592   public:
    00593
    00594     typedef UInt32 argument_type; ///< STL required functor argument type
    00595     typedef UInt32 result_type; ///< STL required functor result type
    00596
    00597         /** Create functor for uniform random integers in the range [0, 2<sup>32</sup>)
    00598             using the given engine.
    00599
    00600             That is, the generated numbers satisfy 0 &lt;= i &lt; 2<sup>32</sup>.
    00601         */
    00602     explicit UniformIntRandomFunctor(Engine & generator = Engine::global() )
    00603     : lower_(0), difference_(0xffffffff), factor_(1),
    00604       generator_(generator),
    00605       useLowBits_(true)
    00606     {}
    00607
    00608         /** Create functor for uniform random integers in the range [<tt>lower</tt>, <tt>upper</tt>]
    00609             using the given engine.
    00610
    00611             That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
    00612             \a useLowBits should be set to <tt>false</tt> when the engine generates
    00613             random numbers whose low bits are significantly less random than the high
    00614             bits. This does not apply to <tt>RandomTT800</tt> and <tt>RandomMT19937</tt>,
    00615             but is necessary for simpler linear congruential generators.
    00616         */
    00617     UniformIntRandomFunctor(UInt32 lower, UInt32 upper,
    00618                             Engine & generator = Engine::global(),
    00619                             bool useLowBits = true)
    00620     : lower_(lower), difference_(upper-lower),
    00621       factor_(Engine::factorForUniformInt(difference_ + 1)),
    00622       generator_(generator),
    00623       useLowBits_(useLowBits)
    00624     {
    00625         vigra_precondition(lower < upper,
    00626           "UniformIntRandomFunctor(): lower bound must be smaller than upper bound.");
    00627     }
    00628
    00629         /** Return a random number as specified in the constructor.
    00630         */
    00631     UInt32 operator()() const
    00632     {
    00633         if(difference_ == 0xffffffff) // lower_ is necessarily 0
    00634             return generator_();
    00635         else if(useLowBits_)
    00636             return generator_.uniformInt(difference_+1) + lower_;
    00637         else
    00638         {
    00639             UInt32 res = generator_() / factor_;
    00640
    00641             // Use rejection method to avoid quantization bias.
    00642             // On average, we will need two raw random numbers to generate one.
    00643             while(res > difference_)
    00644                 res = generator_() / factor_;
    00645             return res + lower_;
    00646         }
    00647     }
    00648
    00649         /** Return a uniformly distributed integer random number in the range [0, <tt>beyond</tt>).
    00650
    00651             That is, 0 &lt;= i &lt; <tt>beyond</tt>. This is a required interface for
    00652             <tt>std::random_shuffle</tt>. It ignores the limits specified
    00653             in the constructor and the flag <tt>useLowBits</tt>.
    00654         */
    00655     UInt32 operator()(UInt32 beyond) const
    00656     {
    00657         if(beyond < 2)
    00658             return 0;
    00659
    00660         return generator_.uniformInt(beyond);
    00661     }
    00662 };
    00663
    00664 template <class Engine>
    00665 class FunctorTraits<UniformIntRandomFunctor<Engine> >
    00666 {
    00667   public:
    00668     typedef UniformIntRandomFunctor<Engine> type;
    00669
    00670     typedef VigraTrueType  isInitializer;
    00671
    00672     typedef VigraTrueType  isUnaryFunctor;
    00673     typedef VigraFalseType isBinaryFunctor;
    00674     typedef VigraFalseType isTernaryFunctor;
    00675
    00676     typedef VigraFalseType isUnaryAnalyser;
    00677     typedef VigraFalseType isBinaryAnalyser;
    00678     typedef VigraFalseType isTernaryAnalyser;
    00679 };
    00680
    00681 /** Functor to create uniformely distributed double-precision random numbers.
    00682
    00683     This functor encapsulates the function <tt>uniform()</tt> of the given random number
    00684     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
    00685     in an STL-compatible interface.
    00686
    00687     <b>Traits defined:</b>
    00688
    00689     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
    00690 */
    00691 template <class Engine = RandomTT800>
    00692 class UniformRandomFunctor
    00693 {
    00694     double offset_, scale_;
    00695     Engine & generator_;
    00696
    00697   public:
    00698
    00699     typedef double result_type; ///< STL required functor result type
    00700
    00701         /** Create functor for uniform random double-precision numbers in the range [0.0, 1.0]
    00702             using the given engine.
    00703
    00704             That is, the generated numbers satisfy 0.0 &lt;= i &lt;= 1.0.
    00705         */
    00706     UniformRandomFunctor(Engine & generator = Engine::global())
    00707     : offset_(0.0),
    00708       scale_(1.0),
    00709       generator_(generator)
    00710     {}
    00711
    00712         /** Create functor for uniform random double-precision numbers in the range [<tt>lower</tt>, <tt>upper</tt>]
    00713             using the given engine.
    00714
    00715             That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
    00716         */
    00717     UniformRandomFunctor(double lower, double upper,
    00718                          Engine & generator = Engine::global())
    00719     : offset_(lower),
    00720       scale_(upper - lower),
    00721       generator_(generator)
    00722     {
    00723         vigra_precondition(lower < upper,
    00724           "UniformRandomFunctor(): lower bound must be smaller than upper bound.");
    00725     }
    00726
    00727         /** Return a random number as specified in the constructor.
    00728         */
    00729     double operator()() const
    00730     {
    00731         return generator_.uniform() * scale_ + offset_;
    00732     }
    00733 };
    00734
    00735 template <class Engine>
    00736 class FunctorTraits<UniformRandomFunctor<Engine> >
    00737 {
    00738   public:
    00739     typedef UniformRandomFunctor<Engine> type;
    00740
    00741     typedef VigraTrueType  isInitializer;
    00742
    00743     typedef VigraFalseType isUnaryFunctor;
    00744     typedef VigraFalseType isBinaryFunctor;
    00745     typedef VigraFalseType isTernaryFunctor;
    00746
    00747     typedef VigraFalseType isUnaryAnalyser;
    00748     typedef VigraFalseType isBinaryAnalyser;
    00749     typedef VigraFalseType isTernaryAnalyser;
    00750 };
    00751
    00752 /** Functor to create normal variate random numbers.
    00753
    00754     This functor encapsulates the function <tt>normal()</tt> of the given random number
    00755     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
    00756     in an STL-compatible interface.
    00757
    00758     <b>Traits defined:</b>
    00759
    00760     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
    00761 */
    00762 template <class Engine = RandomTT800>
    00763 struct NormalRandomFunctor
    00764 {
    00765     double mean_, stddev_;
    00766     Engine & generator_;
    00767
    00768   public:
    00769
    00770     typedef double result_type; ///< STL required functor result type
    00771
    00772         /** Create functor for standard normal random numbers
    00773             using the given engine.
    00774
    00775             That is, mean is 0.0 and standard deviation is 1.0.
    00776         */
    00777     NormalRandomFunctor(Engine & generator = Engine::global())
    00778     : mean_(0.0),
    00779       stddev_(1.0),
    00780       generator_(generator)
    00781     {}
    00782
    00783         /** Create functor for normal random numbers with goven mean and standard deviation
    00784             using the given engine.
    00785         */
    00786     NormalRandomFunctor(double mean, double stddev,
    00787                         Engine & generator = Engine::global())
    00788     : mean_(mean),
    00789       stddev_(stddev),
    00790       generator_(generator)
    00791     {
    00792         vigra_precondition(stddev > 0.0,
    00793           "NormalRandomFunctor(): standard deviation must be positive.");
    00794     }
    00795
    00796         /** Return a random number as specified in the constructor.
    00797         */
    00798     double operator()() const
    00799     {
    00800         return generator_.normal() * stddev_ + mean_;
    00801     }
    00802
    00803 };
    00804
    00805 template <class Engine>
    00806 class FunctorTraits<NormalRandomFunctor<Engine> >
    00807 {
    00808   public:
    00809     typedef UniformRandomFunctor<Engine>  type;
    00810
    00811     typedef VigraTrueType  isInitializer;
    00812
    00813     typedef VigraFalseType isUnaryFunctor;
    00814     typedef VigraFalseType isBinaryFunctor;
    00815     typedef VigraFalseType isTernaryFunctor;
    00816
    00817     typedef VigraFalseType isUnaryAnalyser;
    00818     typedef VigraFalseType isBinaryAnalyser;
    00819     typedef VigraFalseType isTernaryAnalyser;
    00820 };
    00821
    00822 //@}
    00823
    00824 } // namespace vigra
    00825
    00826 #endif // VIGRA_RANDOM_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python VIGRA 1.6.0 (13 Aug 2008)