[ 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://kogs-www.informatik.uni-hamburg.de/~koethe/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();
00144             
00145         UInt32 y = state_[current_++];
00146         y ^= (y << 7) & 0x2b5b2500; 
00147         y ^= (y << 15) & 0xdb8b0000; 
00148         return y ^ (y >> 16);
00149     }
00150     
00151     void generateNumbers() const;
00152 
00153     void seedImpl(RandomSeedTag)
00154     {
00155         seed(RandomSeed, *this);
00156     }
00157 
00158     void seedImpl(UInt32 theSeed)
00159     {
00160         seed(theSeed, *this);
00161     }
00162     
00163     template<class Iterator>
00164     void seedImpl(Iterator init, UInt32 length)
00165     {
00166         seed(init, length, *this);
00167     }
00168 };
00169 
00170 void RandomState<TT800>::generateNumbers() const
00171 {
00172     UInt32 mag01[2]= { 0x0, 0x8ebfd028 };
00173 
00174     for(UInt32 i=0; i<N-M; ++i)
00175     {
00176         state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
00177     }
00178     for (UInt32 i=N-M; i<N; ++i) 
00179     {
00180         state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
00181     }
00182     current_ = 0;
00183 }
00184 
00185     /* Mersenne twister MT19937 by M. Matsumoto */
00186 template<>
00187 struct RandomState<MT19937>
00188 {
00189     static const UInt32 N = 624, M = 397;
00190     
00191     mutable UInt32 state_[N];
00192     mutable UInt32 current_;
00193                    
00194     RandomState()
00195     : current_(0)
00196     {
00197         seed(19650218U, *this);
00198     }
00199 
00200   protected:  
00201 
00202     UInt32 get() const
00203     {
00204         if(current_ == N)
00205             generateNumbers();
00206             
00207         UInt32 x = state_[current_++];
00208         x ^= (x >> 11);
00209         x ^= (x << 7) & 0x9D2C5680U;
00210         x ^= (x << 15) & 0xEFC60000U;
00211         return x ^ (x >> 18);
00212     }
00213     
00214     void generateNumbers() const;
00215 
00216     static UInt32 twiddle(UInt32 u, UInt32 v) 
00217     {
00218         return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
00219                 ^ ((v & 1U) ? 0x9908B0DFU : 0x0U);
00220     }
00221 
00222     void seedImpl(RandomSeedTag)
00223     {
00224         seed(RandomSeed, *this);
00225         generateNumbers();
00226     }
00227 
00228     void seedImpl(UInt32 theSeed)
00229     {
00230         seed(theSeed, *this);
00231         generateNumbers();
00232     }
00233     
00234     template<class Iterator>
00235     void seedImpl(Iterator init, UInt32 length)
00236     {
00237         seed(19650218U, *this);
00238         seed(init, length, *this);
00239         generateNumbers();
00240     }
00241 };
00242 
00243 void RandomState<MT19937>::generateNumbers() const
00244 {
00245     for (int i = 0; i < (N - M); ++i)
00246     {
00247         state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]);
00248     }
00249     for (int i = N - M; i < (N - 1); ++i)
00250     {
00251         state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]);
00252     }
00253     state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]);
00254     current_ = 0;
00255 }
00256 
00257 } // namespace detail
00258 
00259 
00260 /** \addtogroup RandomNumberGeneration Random Number Generation
00261 
00262      High-quality random number generators and related functors.
00263 */
00264 //@{
00265 
00266 /** Generic random number generator.
00267 
00268     The actual generator is passed in the template argument <tt>Engine</tt>. Two generators
00269     are currently available:
00270     <ul>
00271     <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.
00272     <li> <tt>RandomTT800</tt>: (default) The Tempered Twister, a simpler predecessor of the Mersenne Twister with period length 2<sup>800</sup>.
00273     </ul>
00274     
00275     Both generators have been designed by <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/eindex.html">Makoto Matsumoto</a>. 
00276     
00277     <b>Traits defined:</b>
00278     
00279     <tt>FunctorTraits<RandomNumberGenerator<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
00280 */
00281 template <class Engine = detail::RandomState<detail::TT800> >
00282 class RandomNumberGenerator
00283 : public Engine
00284 {
00285     mutable double normalCached_;
00286     mutable bool normalCachedValid_;
00287     
00288   public:
00289   
00290         /** Create a new random generator object with standard seed.
00291             
00292             Due to standard seeding, the random numbers generated will always be the same. 
00293             This is useful for debugging.
00294         */
00295     RandomNumberGenerator()
00296     : normalCached_(0.0),
00297       normalCachedValid_(false)
00298     {}
00299   
00300         /** Create a new random generator object with a random seed.
00301         
00302             The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
00303             values.
00304         
00305             <b>Usage:</b>
00306             \code
00307             RandomNumberGenerator<> rnd = RandomNumberGenerator<>(RandomSeed);
00308             \endcode
00309         */
00310     RandomNumberGenerator(RandomSeedTag)
00311     : normalCached_(0.0),
00312       normalCachedValid_(false)
00313     {
00314         this->seedImpl(RandomSeed);
00315     }
00316   
00317         /** Create a new random generator object from the given seed.
00318             
00319             The same seed will always produce identical random sequences.
00320         */
00321     RandomNumberGenerator(UInt32 theSeed)
00322     : normalCached_(0.0),
00323       normalCachedValid_(false)
00324     {
00325         this->seedImpl(theSeed);
00326     }
00327 
00328         /** Create a new random generator object from the given seed sequence.
00329             
00330             Longer seed sequences lead to better initialization in the sense that the generator's 
00331             state space is covered much better than is possible with 32-bit seeds alone.
00332         */
00333     template<class Iterator>
00334     RandomNumberGenerator(Iterator init, UInt32 length)
00335     : normalCached_(0.0),
00336       normalCachedValid_(false)
00337     {
00338         this->seedImpl(init, length);
00339     }
00340 
00341   
00342         /** Re-initialize the random generator object with a random seed.
00343         
00344             The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
00345             values.
00346         
00347             <b>Usage:</b>
00348             \code
00349             RandomNumberGenerator<> rnd = ...;
00350             ...
00351             rnd.seed(RandomSeed);
00352             \endcode
00353         */
00354     void seed(RandomSeedTag)
00355     {
00356         this->seedImpl(RandomSeed);
00357         normalCachedValid_ = false;
00358     }
00359 
00360         /** Re-initialize the random generator object from the given seed.
00361             
00362             The same seed will always produce identical random sequences.
00363         */
00364     void seed(UInt32 theSeed)
00365     {
00366         this->seedImpl(theSeed);
00367         normalCachedValid_ = false;
00368     }
00369 
00370         /** Re-initialize the random generator object from the given seed sequence.
00371             
00372             Longer seed sequences lead to better initialization in the sense that the generator's 
00373             state space is covered much better than is possible with 32-bit seeds alone.
00374         */
00375     template<class Iterator>
00376     void seed(Iterator init, UInt32 length)
00377     {
00378         this->seedImpl(init, length);
00379         normalCachedValid_ = false;
00380     }
00381 
00382         /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
00383             
00384             That is, 0 &lt;= i &lt; 2<sup>32</sup>. 
00385         */
00386     UInt32 operator()() const
00387     {
00388         return this->get();
00389     }
00390 
00391         /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
00392             
00393             That is, 0 &lt;= i &lt; 2<sup>32</sup>. 
00394         */
00395     UInt32 uniformInt() const
00396     {
00397         return this->get();
00398     }
00399 
00400 
00401 #if 0 // difficult implementation necessary if low bits are not sufficiently random
00402         // in [0,beyond)
00403     UInt32 uniformInt(UInt32 beyond) const
00404     {
00405         if(beyond < 2)
00406             return 0;
00407 
00408         UInt32 factor = factorForUniformInt(beyond);
00409         UInt32 res = this->get() / factor;
00410 
00411         // Use rejection method to avoid quantization bias.
00412         // On average, we will need two raw random numbers to generate one.
00413         while(res >= beyond)
00414             res = this->get() / factor;
00415         return res;
00416     }
00417 #endif /* #if 0 */
00418 
00419         /** Return a uniformly distributed integer random number in [0, <tt>beyond</tt>).
00420             
00421             That is, 0 &lt;= i &lt; <tt>beyond</tt>. 
00422         */
00423     UInt32 uniformInt(UInt32 beyond) const
00424     {
00425         // in [0,beyond) -- simple implementation when low bits are sufficiently random, which is
00426         // the case for TT800 and MT19937
00427         if(beyond < 2)
00428             return 0;
00429 
00430         UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond;
00431         UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder;
00432         UInt32 res = this->get();
00433 
00434         // Use rejection method to avoid quantization bias.
00435         // We will need two raw random numbers in amortized worst case.
00436         while(res > lastSafeValue)
00437             res = this->get();
00438         return res % beyond;
00439     }
00440     
00441         /** Return a uniformly distributed double-precision random number in [0.0, 1.0).
00442             
00443             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 
00444             create this number).
00445         */
00446     double uniform53() const
00447     {
00448         // make full use of the entire 53-bit mantissa of a double, by Isaku Wada
00449         return ( (this->get() >> 5) * 67108864.0 + (this->get() >> 6)) * (1.0/9007199254740992.0); 
00450     }
00451     
00452         /** Return a uniformly distributed double-precision random number in [0.0, 1.0].
00453             
00454             That is, 0.0 &lt;= i &lt;= 1.0. This nuber is computed by <tt>uniformInt()</tt> / 2<sup>32</sup>, 
00455             so it has effectively only 32 random bits. 
00456         */
00457     double uniform() const
00458     {
00459         return (double)this->get() / 4294967295.0;
00460     }
00461 
00462         /** Return a uniformly distributed double-precision random number in [lower, upper].
00463            
00464             That is, <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>. This number is computed 
00465             from <tt>uniform()</tt>, so it has effectively only 32 random bits. 
00466         */
00467     double uniform(double lower, double upper) const
00468     {
00469         vigra_precondition(lower < upper,
00470           "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound."); 
00471         return uniform() * (upper-lower) + lower;
00472     }
00473 
00474         /** Return a standard normal variate (Gaussian) random number.
00475            
00476             Mean is zero, standard deviation is 1.0. It uses the polar form of the 
00477             Box-Muller transform.
00478         */
00479     double normal() const;
00480     
00481         /** Return a normal variate (Gaussian) random number with the given mean and standard deviation.
00482            
00483             It uses the polar form of the Box-Muller transform.
00484         */
00485     double normal(double mean, double stddev) const
00486     {
00487         vigra_precondition(stddev > 0.0,
00488           "RandomNumberGenerator::normal(): standard deviation must be positive."); 
00489         return normal()*stddev + mean;
00490     }
00491     
00492         /** Access the global (program-wide) instance of the present random number generator.
00493         
00494             Normally, you will create a local generator by one of the constructor calls. But sometimes
00495             it is useful to have all program parts access the same generator.
00496         */
00497     static RandomNumberGenerator & global()
00498     {
00499         static RandomNumberGenerator generator;
00500         return generator;
00501     }
00502 
00503     static UInt32 factorForUniformInt(UInt32 range)
00504     {
00505         return (range > 2147483648U || range == 0)
00506                      ? 1
00507                      : 2*(2147483648U / ceilPower2(range));
00508     }
00509 };
00510 
00511 template <class Engine>
00512 double RandomNumberGenerator<Engine>::normal() const
00513 {
00514     if(normalCachedValid_)
00515     {
00516         normalCachedValid_ = false;
00517         return normalCached_;
00518     }
00519     else
00520     {
00521         double x1, x2, w;
00522         do 
00523         {
00524              x1 = uniform(-1.0, 1.0);
00525              x2 = uniform(-1.0, 1.0);
00526              w = x1 * x1 + x2 * x2;
00527         } 
00528         while ( w > 1.0 || w == 0.0);
00529         
00530         w = std::sqrt( -2.0 * std::log( w )  / w );
00531 
00532         normalCached_ = x2 * w;
00533         normalCachedValid_ = true;
00534 
00535         return x1 * w;
00536     }
00537 }
00538 
00539     /** Shorthand for the TT800 random number generator class.
00540     */
00541 typedef RandomNumberGenerator<>  RandomTT800; 
00542 
00543     /** Shorthand for the MT19937 random number generator class.
00544     */
00545 typedef RandomNumberGenerator<detail::RandomState<detail::MT19937> > RandomMT19937;
00546 
00547     /** Access the global (program-wide) instance of the TT800 random number generator.
00548     */
00549 inline RandomTT800   & randomTT800()   { return RandomTT800::global(); }
00550 
00551     /** Access the global (program-wide) instance of the MT19937 random number generator.
00552     */
00553 inline RandomMT19937 & randomMT19937() { return RandomMT19937::global(); }
00554 
00555 template <class Engine>
00556 class FunctorTraits<RandomNumberGenerator<Engine> >
00557 {
00558   public:
00559     typedef RandomNumberGenerator<Engine> type;
00560     
00561     typedef VigraTrueType  isInitializer;
00562     
00563     typedef VigraFalseType isUnaryFunctor;
00564     typedef VigraFalseType isBinaryFunctor;
00565     typedef VigraFalseType isTernaryFunctor;
00566     
00567     typedef VigraFalseType isUnaryAnalyser;
00568     typedef VigraFalseType isBinaryAnalyser;
00569     typedef VigraFalseType isTernaryAnalyser;
00570 };
00571 
00572 
00573 /** Functor to create uniformely distributed integer random numbers.
00574 
00575     This functor encapsulates the appropriate functions of the given random number
00576     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
00577     in an STL-compatible interface. 
00578     
00579     <b>Traits defined:</b>
00580     
00581     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> and
00582     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isUnaryFunctor</tt> are true (<tt>VigraTrueType</tt>).
00583 */
00584 template <class Engine = RandomTT800>
00585 class UniformIntRandomFunctor
00586 {
00587     UInt32 lower_, difference_, factor_;
00588     Engine & generator_;
00589     bool useLowBits_;
00590 
00591   public:
00592   
00593     typedef UInt32 argument_type; ///< STL required functor argument type
00594     typedef UInt32 result_type; ///< STL required functor result type
00595 
00596         /** Create functor for uniform random integers in the range [0, 2<sup>32</sup>) 
00597             using the given engine.
00598             
00599             That is, the generated numbers satisfy 0 &lt;= i &lt; 2<sup>32</sup>.
00600         */
00601     explicit UniformIntRandomFunctor(Engine & generator = Engine::global() )
00602     : lower_(0), difference_(0xffffffff), factor_(1),
00603       generator_(generator),
00604       useLowBits_(true)
00605     {}
00606     
00607         /** Create functor for uniform random integers in the range [<tt>lower</tt>, <tt>upper</tt>]
00608             using the given engine.
00609             
00610             That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
00611             \a useLowBits should be set to <tt>false</tt> when the engine generates
00612             random numbers whose low bits are significantly less random than the high
00613             bits. This does not apply to <tt>RandomTT800</tt> and <tt>RandomMT19937</tt>,
00614             but is necessary for simpler linear congruential generators.
00615         */
00616     UniformIntRandomFunctor(UInt32 lower, UInt32 upper, 
00617                             Engine & generator = Engine::global(),
00618                             bool useLowBits = true)
00619     : lower_(lower), difference_(upper-lower), 
00620       factor_(Engine::factorForUniformInt(difference_ + 1)),
00621       generator_(generator),
00622       useLowBits_(useLowBits)
00623     {
00624         vigra_precondition(lower < upper,
00625           "UniformIntRandomFunctor(): lower bound must be smaller than upper bound."); 
00626     }
00627     
00628         /** Return a random number as specified in the constructor.
00629         */
00630     UInt32 operator()() const
00631     {
00632         if(difference_ == 0xffffffff) // lower_ is necessarily 0
00633             return generator_();
00634         else if(useLowBits_)
00635             return generator_.uniformInt(difference_+1) + lower_;
00636         else
00637         {
00638             UInt32 res = generator_() / factor_;
00639 
00640             // Use rejection method to avoid quantization bias.
00641             // On average, we will need two raw random numbers to generate one.
00642             while(res > difference_)
00643                 res = generator_() / factor_;
00644             return res + lower_;
00645         }
00646     }
00647 
00648         /** Return a uniformly distributed integer random number in the range [0, <tt>beyond</tt>).
00649         
00650             That is, 0 &lt;= i &lt; <tt>beyond</tt>. This is a required interface for 
00651             <tt>std::random_shuffle</tt>. It ignores the limits specified 
00652             in the constructor and the flag <tt>useLowBits</tt>.
00653         */
00654     UInt32 operator()(UInt32 beyond) const
00655     {
00656         if(beyond < 2)
00657             return 0;
00658 
00659         return generator_.uniformInt(beyond);
00660     }
00661 };
00662 
00663 template <class Engine>
00664 class FunctorTraits<UniformIntRandomFunctor<Engine> >
00665 {
00666   public:
00667     typedef UniformIntRandomFunctor<Engine> type;
00668     
00669     typedef VigraTrueType  isInitializer;
00670     
00671     typedef VigraTrueType  isUnaryFunctor;
00672     typedef VigraFalseType isBinaryFunctor;
00673     typedef VigraFalseType isTernaryFunctor;
00674     
00675     typedef VigraFalseType isUnaryAnalyser;
00676     typedef VigraFalseType isBinaryAnalyser;
00677     typedef VigraFalseType isTernaryAnalyser;
00678 };
00679 
00680 /** Functor to create uniformely distributed double-precision random numbers.
00681 
00682     This functor encapsulates the function <tt>uniform()</tt> of the given random number
00683     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
00684     in an STL-compatible interface. 
00685     
00686     <b>Traits defined:</b>
00687     
00688     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
00689 */
00690 template <class Engine = RandomTT800>
00691 class UniformRandomFunctor
00692 {
00693     double offset_, scale_;
00694     Engine & generator_;
00695 
00696   public:
00697   
00698     typedef double result_type; ///< STL required functor result type
00699 
00700         /** Create functor for uniform random double-precision numbers in the range [0.0, 1.0] 
00701             using the given engine.
00702             
00703             That is, the generated numbers satisfy 0.0 &lt;= i &lt;= 1.0.
00704         */
00705     UniformRandomFunctor(Engine & generator = Engine::global())
00706     : offset_(0.0),
00707       scale_(1.0),
00708       generator_(generator)
00709     {}
00710 
00711         /** Create functor for uniform random double-precision numbers in the range [<tt>lower</tt>, <tt>upper</tt>]
00712             using the given engine.
00713             
00714             That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
00715         */
00716     UniformRandomFunctor(double lower, double upper, 
00717                          Engine & generator = Engine::global())
00718     : offset_(lower),
00719       scale_(upper - lower),
00720       generator_(generator)
00721     {
00722         vigra_precondition(lower < upper,
00723           "UniformRandomFunctor(): lower bound must be smaller than upper bound."); 
00724     }
00725     
00726         /** Return a random number as specified in the constructor.
00727         */
00728     double operator()() const
00729     {
00730         return generator_.uniform() * scale_ + offset_;
00731     }
00732 };
00733 
00734 template <class Engine>
00735 class FunctorTraits<UniformRandomFunctor<Engine> >
00736 {
00737   public:
00738     typedef UniformRandomFunctor<Engine> type;
00739     
00740     typedef VigraTrueType  isInitializer;
00741     
00742     typedef VigraFalseType isUnaryFunctor;
00743     typedef VigraFalseType isBinaryFunctor;
00744     typedef VigraFalseType isTernaryFunctor;
00745     
00746     typedef VigraFalseType isUnaryAnalyser;
00747     typedef VigraFalseType isBinaryAnalyser;
00748     typedef VigraFalseType isTernaryAnalyser;
00749 };
00750 
00751 /** Functor to create normal variate random numbers.
00752 
00753     This functor encapsulates the function <tt>normal()</tt> of the given random number
00754     <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
00755     in an STL-compatible interface. 
00756     
00757     <b>Traits defined:</b>
00758     
00759     <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
00760 */
00761 template <class Engine = RandomTT800>
00762 struct NormalRandomFunctor
00763 {
00764     double mean_, stddev_;
00765     Engine & generator_;
00766 
00767   public:
00768   
00769     typedef double result_type; ///< STL required functor result type
00770 
00771         /** Create functor for standard normal random numbers 
00772             using the given engine.
00773             
00774             That is, mean is 0.0 and standard deviation is 1.0.
00775         */
00776     NormalRandomFunctor(Engine & generator = Engine::global())
00777     : mean_(0.0),
00778       stddev_(1.0),
00779       generator_(generator)
00780     {}
00781 
00782         /** Create functor for normal random numbers with goven mean and standard deviation
00783             using the given engine.
00784         */
00785     NormalRandomFunctor(double mean, double stddev, 
00786                         Engine & generator = Engine::global())
00787     : mean_(mean),
00788       stddev_(stddev),
00789       generator_(generator)
00790     {
00791         vigra_precondition(stddev > 0.0,
00792           "NormalRandomFunctor(): standard deviation must be positive."); 
00793     }
00794     
00795         /** Return a random number as specified in the constructor.
00796         */
00797     double operator()() const
00798     {
00799         return generator_.normal() * stddev_ + mean_;
00800     }
00801 
00802 };
00803 
00804 template <class Engine>
00805 class FunctorTraits<NormalRandomFunctor<Engine> >
00806 {
00807   public:
00808     typedef UniformRandomFunctor<Engine>  type;
00809     
00810     typedef VigraTrueType  isInitializer;
00811     
00812     typedef VigraFalseType isUnaryFunctor;
00813     typedef VigraFalseType isBinaryFunctor;
00814     typedef VigraFalseType isTernaryFunctor;
00815     
00816     typedef VigraFalseType isUnaryAnalyser;
00817     typedef VigraFalseType isBinaryAnalyser;
00818     typedef VigraFalseType isTernaryAnalyser;
00819 };
00820 
00821 //@}
00822 
00823 } // namespace vigra 
00824 
00825 #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 (5 Nov 2009)