[ 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)