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

random.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2008 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 
37 #ifndef VIGRA_RANDOM_HXX
38 #define VIGRA_RANDOM_HXX
39 
40 #include "mathutil.hxx"
41 #include "functortraits.hxx"
42 
43 #include <ctime>
44 
45 namespace vigra {
46 
47 enum RandomSeedTag { RandomSeed };
48 
49 namespace detail {
50 
51 enum RandomEngineTag { TT800, MT19937 };
52 
53 template<RandomEngineTag EngineTag>
54 struct RandomState;
55 
56 template <RandomEngineTag EngineTag>
57 void seed(UInt32 theSeed, RandomState<EngineTag> & engine)
58 {
59  engine.state_[0] = theSeed;
60  for(UInt32 i=1; i<RandomState<EngineTag>::N; ++i)
61  {
62  engine.state_[i] = 1812433253U * (engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) + i;
63  }
64 }
65 
66 template <class Iterator, RandomEngineTag EngineTag>
67 void seed(Iterator init, UInt32 key_length, RandomState<EngineTag> & engine)
68 {
69  const UInt32 N = RandomState<EngineTag>::N;
70  int k = (int)std::max(N, key_length);
71  UInt32 i = 1, j = 0;
72  Iterator data = init;
73  for (; k; --k)
74  {
75  engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1664525U))
76  + *data + j; /* non linear */
77  ++i; ++j; ++data;
78 
79  if (i >= N)
80  {
81  engine.state_[0] = engine.state_[N-1];
82  i=1;
83  }
84  if (j>=key_length)
85  {
86  j=0;
87  data = init;
88  }
89  }
90 
91  for (k=N-1; k; --k)
92  {
93  engine.state_[i] = (engine.state_[i] ^ ((engine.state_[i-1] ^ (engine.state_[i-1] >> 30)) * 1566083941U))
94  - i; /* non linear */
95  ++i;
96  if (i>=N)
97  {
98  engine.state_[0] = engine.state_[N-1];
99  i=1;
100  }
101  }
102 
103  engine.state_[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
104 }
105 
106 template <RandomEngineTag EngineTag>
107 void seed(RandomSeedTag, RandomState<EngineTag> & engine)
108 {
109  static UInt32 globalCount = 0;
110  UInt32 init[3] = { (UInt32)time(0), (UInt32)clock(), ++globalCount };
111  seed(init, 3, engine);
112 }
113 
114 
115  /* Tempered twister TT800 by M. Matsumoto */
116 template<>
117 struct RandomState<TT800>
118 {
119  static const UInt32 N = 25, M = 7;
120 
121  mutable UInt32 state_[N];
122  mutable UInt32 current_;
123 
124  RandomState()
125  : current_(0)
126  {
127  UInt32 seeds[N] = {
128  0x95f24dab, 0x0b685215, 0xe76ccae7, 0xaf3ec239, 0x715fad23,
129  0x24a590ad, 0x69e4b5ef, 0xbf456141, 0x96bc1b7b, 0xa7bdf825,
130  0xc1de75b7, 0x8858a9c9, 0x2da87693, 0xb657f9dd, 0xffdc8a9f,
131  0x8121da71, 0x8b823ecb, 0x885d05f5, 0x4e20cd47, 0x5a9ad5d9,
132  0x512c0c03, 0xea857ccd, 0x4cc1d30f, 0x8891a8a1, 0xa6b7aadb
133  };
134 
135  for(UInt32 i=0; i<N; ++i)
136  state_[i] = seeds[i];
137  }
138 
139  protected:
140 
141  UInt32 get() const
142  {
143  if(current_ == N)
144  generateNumbers<void>();
145 
146  UInt32 y = state_[current_++];
147  y ^= (y << 7) & 0x2b5b2500;
148  y ^= (y << 15) & 0xdb8b0000;
149  return y ^ (y >> 16);
150  }
151 
152  template <class DUMMY>
153  void generateNumbers() const;
154 
155  void seedImpl(RandomSeedTag)
156  {
157  seed(RandomSeed, *this);
158  }
159 
160  void seedImpl(UInt32 theSeed)
161  {
162  seed(theSeed, *this);
163  }
164 
165  template<class Iterator>
166  void seedImpl(Iterator init, UInt32 length)
167  {
168  seed(init, length, *this);
169  }
170 };
171 
172 template <class DUMMY>
173 void RandomState<TT800>::generateNumbers() const
174 {
175  UInt32 mag01[2]= { 0x0, 0x8ebfd028 };
176 
177  for(UInt32 i=0; i<N-M; ++i)
178  {
179  state_[i] = state_[i+M] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
180  }
181  for (UInt32 i=N-M; i<N; ++i)
182  {
183  state_[i] = state_[i+(M-N)] ^ (state_[i] >> 1) ^ mag01[state_[i] % 2];
184  }
185  current_ = 0;
186 }
187 
188  /* Mersenne twister MT19937 by M. Matsumoto */
189 template<>
190 struct RandomState<MT19937>
191 {
192  static const UInt32 N = 624, M = 397;
193 
194  mutable UInt32 state_[N];
195  mutable UInt32 current_;
196 
197  RandomState()
198  : current_(0)
199  {
200  seed(19650218U, *this);
201  }
202 
203  protected:
204 
205  UInt32 get() const
206  {
207  if(current_ == N)
208  generateNumbers<void>();
209 
210  UInt32 x = state_[current_++];
211  x ^= (x >> 11);
212  x ^= (x << 7) & 0x9D2C5680U;
213  x ^= (x << 15) & 0xEFC60000U;
214  return x ^ (x >> 18);
215  }
216 
217  template <class DUMMY>
218  void generateNumbers() const;
219 
220  static UInt32 twiddle(UInt32 u, UInt32 v)
221  {
222  return (((u & 0x80000000U) | (v & 0x7FFFFFFFU)) >> 1)
223  ^ ((v & 1U) ? 0x9908B0DFU : 0x0U);
224  }
225 
226  void seedImpl(RandomSeedTag)
227  {
228  seed(RandomSeed, *this);
229  generateNumbers<void>();
230  }
231 
232  void seedImpl(UInt32 theSeed)
233  {
234  seed(theSeed, *this);
235  generateNumbers<void>();
236  }
237 
238  template<class Iterator>
239  void seedImpl(Iterator init, UInt32 length)
240  {
241  seed(19650218U, *this);
242  seed(init, length, *this);
243  generateNumbers<void>();
244  }
245 };
246 
247 template <class DUMMY>
248 void RandomState<MT19937>::generateNumbers() const
249 {
250  for (unsigned int i = 0; i < (N - M); ++i)
251  {
252  state_[i] = state_[i + M] ^ twiddle(state_[i], state_[i + 1]);
253  }
254  for (unsigned int i = N - M; i < (N - 1); ++i)
255  {
256  state_[i] = state_[i + M - N] ^ twiddle(state_[i], state_[i + 1]);
257  }
258  state_[N - 1] = state_[M - 1] ^ twiddle(state_[N - 1], state_[0]);
259  current_ = 0;
260 }
261 
262 } // namespace detail
263 
264 
265 /** \addtogroup RandomNumberGeneration Random Number Generation
266 
267  High-quality random number generators and related functors.
268 */
269 //@{
270 
271 /** Generic random number generator.
272 
273  The actual generator is passed in the template argument <tt>Engine</tt>. Two generators
274  are currently available:
275  <ul>
276  <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.
277  <li> <tt>RandomTT800</tt>: (default) The Tempered Twister, a simpler predecessor of the Mersenne Twister with period length 2<sup>800</sup>.
278  </ul>
279 
280  Both generators have been designed by <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/eindex.html">Makoto Matsumoto</a>.
281 
282  <b>Traits defined:</b>
283 
284  <tt>FunctorTraits<RandomNumberGenerator<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
285 */
286 template <class Engine = detail::RandomState<detail::TT800> >
288 : public Engine
289 {
290  mutable double normalCached_;
291  mutable bool normalCachedValid_;
292 
293  public:
294 
295  /** Create a new random generator object with standard seed.
296 
297  Due to standard seeding, the random numbers generated will always be the same.
298  This is useful for debugging.
299  */
301  : normalCached_(0.0),
302  normalCachedValid_(false)
303  {}
304 
305  /** Create a new random generator object with a random seed.
306 
307  The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
308  values.
309 
310  <b>Usage:</b>
311  \code
312  RandomNumberGenerator<> rnd = RandomNumberGenerator<>(RandomSeed);
313  \endcode
314  */
315  RandomNumberGenerator(RandomSeedTag)
316  : normalCached_(0.0),
317  normalCachedValid_(false)
318  {
319  this->seedImpl(RandomSeed);
320  }
321 
322  /** Create a new random generator object from the given seed.
323 
324  The same seed will always produce identical random sequences.
325  */
327  : normalCached_(0.0),
328  normalCachedValid_(false)
329  {
330  this->seedImpl(theSeed);
331  }
332 
333  /** Create a new random generator object from the given seed sequence.
334 
335  Longer seed sequences lead to better initialization in the sense that the generator's
336  state space is covered much better than is possible with 32-bit seeds alone.
337  */
338  template<class Iterator>
339  RandomNumberGenerator(Iterator init, UInt32 length)
340  : normalCached_(0.0),
341  normalCachedValid_(false)
342  {
343  this->seedImpl(init, length);
344  }
345 
346 
347  /** Re-initialize the random generator object with a random seed.
348 
349  The seed is obtained from the machines current <tt>clock()</tt> and <tt>time()</tt>
350  values.
351 
352  <b>Usage:</b>
353  \code
354  RandomNumberGenerator<> rnd = ...;
355  ...
356  rnd.seed(RandomSeed);
357  \endcode
358  */
359  void seed(RandomSeedTag)
360  {
361  this->seedImpl(RandomSeed);
362  normalCachedValid_ = false;
363  }
364 
365  /** Re-initialize the random generator object from the given seed.
366 
367  The same seed will always produce identical random sequences.
368  */
369  void seed(UInt32 theSeed)
370  {
371  this->seedImpl(theSeed);
372  normalCachedValid_ = false;
373  }
374 
375  /** Re-initialize the random generator object from the given seed sequence.
376 
377  Longer seed sequences lead to better initialization in the sense that the generator's
378  state space is covered much better than is possible with 32-bit seeds alone.
379  */
380  template<class Iterator>
381  void seed(Iterator init, UInt32 length)
382  {
383  this->seedImpl(init, length);
384  normalCachedValid_ = false;
385  }
386 
387  /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
388 
389  That is, 0 &lt;= i &lt; 2<sup>32</sup>.
390  */
392  {
393  return this->get();
394  }
395 
396  /** Return a uniformly distributed integer random number in [0, 2<sup>32</sup>).
397 
398  That is, 0 &lt;= i &lt; 2<sup>32</sup>.
399  */
401  {
402  return this->get();
403  }
404 
405 
406 #if 0 // difficult implementation necessary if low bits are not sufficiently random
407  // in [0,beyond)
408  UInt32 uniformInt(UInt32 beyond) const
409  {
410  if(beyond < 2)
411  return 0;
412 
413  UInt32 factor = factorForUniformInt(beyond);
414  UInt32 res = this->get() / factor;
415 
416  // Use rejection method to avoid quantization bias.
417  // On average, we will need two raw random numbers to generate one.
418  while(res >= beyond)
419  res = this->get() / factor;
420  return res;
421  }
422 #endif /* #if 0 */
423 
424  /** Return a uniformly distributed integer random number in [0, <tt>beyond</tt>).
425 
426  That is, 0 &lt;= i &lt; <tt>beyond</tt>.
427  */
428  UInt32 uniformInt(UInt32 beyond) const
429  {
430  // in [0,beyond) -- simple implementation when low bits are sufficiently random, which is
431  // the case for TT800 and MT19937
432  if(beyond < 2)
433  return 0;
434 
435  UInt32 remainder = (NumericTraits<UInt32>::max() - beyond + 1) % beyond;
436  UInt32 lastSafeValue = NumericTraits<UInt32>::max() - remainder;
437  UInt32 res = this->get();
438 
439  // Use rejection method to avoid quantization bias.
440  // We will need two raw random numbers in amortized worst case.
441  while(res > lastSafeValue)
442  res = this->get();
443  return res % beyond;
444  }
445 
446  /** Return a uniformly distributed double-precision random number in [0.0, 1.0).
447 
448  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
449  create this number).
450  */
451  double uniform53() const
452  {
453  // make full use of the entire 53-bit mantissa of a double, by Isaku Wada
454  return ( (this->get() >> 5) * 67108864.0 + (this->get() >> 6)) * (1.0/9007199254740992.0);
455  }
456 
457  /** Return a uniformly distributed double-precision random number in [0.0, 1.0].
458 
459  That is, 0.0 &lt;= i &lt;= 1.0. This nuber is computed by <tt>uniformInt()</tt> / 2<sup>32</sup>,
460  so it has effectively only 32 random bits.
461  */
462  double uniform() const
463  {
464  return (double)this->get() / 4294967295.0;
465  }
466 
467  /** Return a uniformly distributed double-precision random number in [lower, upper].
468 
469  That is, <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>. This number is computed
470  from <tt>uniform()</tt>, so it has effectively only 32 random bits.
471  */
472  double uniform(double lower, double upper) const
473  {
474  vigra_precondition(lower < upper,
475  "RandomNumberGenerator::uniform(): lower bound must be smaller than upper bound.");
476  return uniform() * (upper-lower) + lower;
477  }
478 
479  /** Return a standard normal variate (Gaussian) random number.
480 
481  Mean is zero, standard deviation is 1.0. It uses the polar form of the
482  Box-Muller transform.
483  */
484  double normal() const;
485 
486  /** Return a normal variate (Gaussian) random number with the given mean and standard deviation.
487 
488  It uses the polar form of the Box-Muller transform.
489  */
490  double normal(double mean, double stddev) const
491  {
492  vigra_precondition(stddev > 0.0,
493  "RandomNumberGenerator::normal(): standard deviation must be positive.");
494  return normal()*stddev + mean;
495  }
496 
497  /** Access the global (program-wide) instance of the present random number generator.
498 
499  Normally, you will create a local generator by one of the constructor calls. But sometimes
500  it is useful to have all program parts access the same generator.
501  */
503  {
504  static RandomNumberGenerator generator;
505  return generator;
506  }
507 
508  static UInt32 factorForUniformInt(UInt32 range)
509  {
510  return (range > 2147483648U || range == 0)
511  ? 1
512  : 2*(2147483648U / ceilPower2(range));
513  }
514 };
515 
516 template <class Engine>
518 {
519  if(normalCachedValid_)
520  {
521  normalCachedValid_ = false;
522  return normalCached_;
523  }
524  else
525  {
526  double x1, x2, w;
527  do
528  {
529  x1 = uniform(-1.0, 1.0);
530  x2 = uniform(-1.0, 1.0);
531  w = x1 * x1 + x2 * x2;
532  }
533  while ( w > 1.0 || w == 0.0);
534 
535  w = std::sqrt( -2.0 * std::log( w ) / w );
536 
537  normalCached_ = x2 * w;
538  normalCachedValid_ = true;
539 
540  return x1 * w;
541  }
542 }
543 
544  /** Shorthand for the TT800 random number generator class.
545  */
547 
548  /** Shorthand for the TT800 random number generator class (same as RandomTT800).
549  */
551 
552  /** Shorthand for the MT19937 random number generator class.
553  */
555 
556  /** Shorthand for the MT19937 random number generator class (same as RandomMT19937).
557  */
559 
560  /** Access the global (program-wide) instance of the TT800 random number generator.
561  */
563 
564  /** Access the global (program-wide) instance of the MT19937 random number generator.
565  */
567 
568 template <class Engine>
569 class FunctorTraits<RandomNumberGenerator<Engine> >
570 {
571  public:
572  typedef RandomNumberGenerator<Engine> type;
573 
574  typedef VigraTrueType isInitializer;
575 
576  typedef VigraFalseType isUnaryFunctor;
577  typedef VigraFalseType isBinaryFunctor;
578  typedef VigraFalseType isTernaryFunctor;
579 
580  typedef VigraFalseType isUnaryAnalyser;
581  typedef VigraFalseType isBinaryAnalyser;
582  typedef VigraFalseType isTernaryAnalyser;
583 };
584 
585 
586 /** Functor to create uniformely distributed integer random numbers.
587 
588  This functor encapsulates the appropriate functions of the given random number
589  <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
590  in an STL-compatible interface.
591 
592  <b>Traits defined:</b>
593 
594  <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> and
595  <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isUnaryFunctor</tt> are true (<tt>VigraTrueType</tt>).
596 */
597 template <class Engine = RandomTT800>
599 {
600  UInt32 lower_, difference_, factor_;
601  Engine const & generator_;
602  bool useLowBits_;
603 
604  public:
605 
606  typedef UInt32 argument_type; ///< STL required functor argument type
607  typedef UInt32 result_type; ///< STL required functor result type
608 
609  /** Create functor for uniform random integers in the range [0, 2<sup>32</sup>)
610  using the given engine.
611 
612  That is, the generated numbers satisfy 0 &lt;= i &lt; 2<sup>32</sup>.
613  */
614  explicit UniformIntRandomFunctor(Engine const & generator = Engine::global() )
615  : lower_(0), difference_(0xffffffff), factor_(1),
616  generator_(generator),
617  useLowBits_(true)
618  {}
619 
620  /** Create functor for uniform random integers in the range [<tt>lower</tt>, <tt>upper</tt>]
621  using the given engine.
622 
623  That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
624  \a useLowBits should be set to <tt>false</tt> when the engine generates
625  random numbers whose low bits are significantly less random than the high
626  bits. This does not apply to <tt>RandomTT800</tt> and <tt>RandomMT19937</tt>,
627  but is necessary for simpler linear congruential generators.
628  */
630  Engine const & generator = Engine::global(),
631  bool useLowBits = true)
632  : lower_(lower), difference_(upper-lower),
633  factor_(Engine::factorForUniformInt(difference_ + 1)),
634  generator_(generator),
635  useLowBits_(useLowBits)
636  {
637  vigra_precondition(lower < upper,
638  "UniformIntRandomFunctor(): lower bound must be smaller than upper bound.");
639  }
640 
641  /** Return a random number as specified in the constructor.
642  */
644  {
645  if(difference_ == 0xffffffff) // lower_ is necessarily 0
646  return generator_();
647  else if(useLowBits_)
648  return generator_.uniformInt(difference_+1) + lower_;
649  else
650  {
651  UInt32 res = generator_() / factor_;
652 
653  // Use rejection method to avoid quantization bias.
654  // On average, we will need two raw random numbers to generate one.
655  while(res > difference_)
656  res = generator_() / factor_;
657  return res + lower_;
658  }
659  }
660 
661  /** Return a uniformly distributed integer random number in the range [0, <tt>beyond</tt>).
662 
663  That is, 0 &lt;= i &lt; <tt>beyond</tt>. This is a required interface for
664  <tt>std::random_shuffle</tt>. It ignores the limits specified
665  in the constructor and the flag <tt>useLowBits</tt>.
666  */
667  UInt32 operator()(UInt32 beyond) const
668  {
669  if(beyond < 2)
670  return 0;
671 
672  return generator_.uniformInt(beyond);
673  }
674 };
675 
676 template <class Engine>
677 class FunctorTraits<UniformIntRandomFunctor<Engine> >
678 {
679  public:
680  typedef UniformIntRandomFunctor<Engine> type;
681 
682  typedef VigraTrueType isInitializer;
683 
684  typedef VigraTrueType isUnaryFunctor;
685  typedef VigraFalseType isBinaryFunctor;
686  typedef VigraFalseType isTernaryFunctor;
687 
688  typedef VigraFalseType isUnaryAnalyser;
689  typedef VigraFalseType isBinaryAnalyser;
690  typedef VigraFalseType isTernaryAnalyser;
691 };
692 
693 /** Functor to create uniformely distributed double-precision random numbers.
694 
695  This functor encapsulates the function <tt>uniform()</tt> of the given random number
696  <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
697  in an STL-compatible interface.
698 
699  <b>Traits defined:</b>
700 
701  <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
702 */
703 template <class Engine = RandomTT800>
705 {
706  double offset_, scale_;
707  Engine const & generator_;
708 
709  public:
710 
711  typedef double result_type; ///< STL required functor result type
712 
713  /** Create functor for uniform random double-precision numbers in the range [0.0, 1.0]
714  using the given engine.
715 
716  That is, the generated numbers satisfy 0.0 &lt;= i &lt;= 1.0.
717  */
718  UniformRandomFunctor(Engine const & generator = Engine::global())
719  : offset_(0.0),
720  scale_(1.0),
721  generator_(generator)
722  {}
723 
724  /** Create functor for uniform random double-precision numbers in the range [<tt>lower</tt>, <tt>upper</tt>]
725  using the given engine.
726 
727  That is, the generated numbers satisfy <tt>lower</tt> &lt;= i &lt;= <tt>upper</tt>.
728  */
729  UniformRandomFunctor(double lower, double upper,
730  Engine & generator = Engine::global())
731  : offset_(lower),
732  scale_(upper - lower),
733  generator_(generator)
734  {
735  vigra_precondition(lower < upper,
736  "UniformRandomFunctor(): lower bound must be smaller than upper bound.");
737  }
738 
739  /** Return a random number as specified in the constructor.
740  */
741  double operator()() const
742  {
743  return generator_.uniform() * scale_ + offset_;
744  }
745 };
746 
747 template <class Engine>
748 class FunctorTraits<UniformRandomFunctor<Engine> >
749 {
750  public:
751  typedef UniformRandomFunctor<Engine> type;
752 
753  typedef VigraTrueType isInitializer;
754 
755  typedef VigraFalseType isUnaryFunctor;
756  typedef VigraFalseType isBinaryFunctor;
757  typedef VigraFalseType isTernaryFunctor;
758 
759  typedef VigraFalseType isUnaryAnalyser;
760  typedef VigraFalseType isBinaryAnalyser;
761  typedef VigraFalseType isTernaryAnalyser;
762 };
763 
764 /** Functor to create normal variate random numbers.
765 
766  This functor encapsulates the function <tt>normal()</tt> of the given random number
767  <tt>Engine</tt> (usually <tt>RandomTT800</tt> or <tt>RandomMT19937</tt>)
768  in an STL-compatible interface.
769 
770  <b>Traits defined:</b>
771 
772  <tt>FunctorTraits<UniformIntRandomFunctor<Engine> >::isInitializer</tt> is true (<tt>VigraTrueType</tt>).
773 */
774 template <class Engine = RandomTT800>
776 {
777  double mean_, stddev_;
778  Engine const & generator_;
779 
780  public:
781 
782  typedef double result_type; ///< STL required functor result type
783 
784  /** Create functor for standard normal random numbers
785  using the given engine.
786 
787  That is, mean is 0.0 and standard deviation is 1.0.
788  */
789  NormalRandomFunctor(Engine const & generator = Engine::global())
790  : mean_(0.0),
791  stddev_(1.0),
792  generator_(generator)
793  {}
794 
795  /** Create functor for normal random numbers with goven mean and standard deviation
796  using the given engine.
797  */
798  NormalRandomFunctor(double mean, double stddev,
799  Engine & generator = Engine::global())
800  : mean_(mean),
801  stddev_(stddev),
802  generator_(generator)
803  {
804  vigra_precondition(stddev > 0.0,
805  "NormalRandomFunctor(): standard deviation must be positive.");
806  }
807 
808  /** Return a random number as specified in the constructor.
809  */
810  double operator()() const
811  {
812  return generator_.normal() * stddev_ + mean_;
813  }
814 
815 };
816 
817 template <class Engine>
818 class FunctorTraits<NormalRandomFunctor<Engine> >
819 {
820  public:
821  typedef UniformRandomFunctor<Engine> type;
822 
823  typedef VigraTrueType isInitializer;
824 
825  typedef VigraFalseType isUnaryFunctor;
826  typedef VigraFalseType isBinaryFunctor;
827  typedef VigraFalseType isTernaryFunctor;
828 
829  typedef VigraFalseType isUnaryAnalyser;
830  typedef VigraFalseType isBinaryAnalyser;
831  typedef VigraFalseType isTernaryAnalyser;
832 };
833 
834 //@}
835 
836 } // namespace vigra
837 
838 #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.1 (Thu Jun 14 2012)