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