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