git.maemo.org Git - opencv/blob - 3rdparty/include/OpenEXR/ImathRandom.h

   1 ///////////////////////////////////////////////////////////////////////////
   2 //
   3 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
   4 // Digital Ltd. LLC
   5 //
   6 // All rights reserved.
   7 //
   8 // Redistribution and use in source and binary forms, with or without
   9 // modification, are permitted provided that the following conditions are
  10 // met:
  11 // *       Redistributions of source code must retain the above copyright
  12 // notice, this list of conditions and the following disclaimer.
  13 // *       Redistributions in binary form must reproduce the above
  14 // copyright notice, this list of conditions and the following disclaimer
  15 // in the documentation and/or other materials provided with the
  16 // distribution.
  17 // *       Neither the name of Industrial Light & Magic nor the names of
  18 // its contributors may be used to endorse or promote products derived
  19 // from this software without specific prior written permission.
  20 //
  21 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  22 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  23 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  24 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  25 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  26 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  27 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  28 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  29 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  30 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  31 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  32 //
  33 ///////////////////////////////////////////////////////////////////////////
  34
  35
  36 #ifndef INCLUDED_IMATHRANDOM_H
  37 #define INCLUDED_IMATHRANDOM_H
  38
  39 //-----------------------------------------------------------------------------
  40 //
  41 //      Generators for uniformly distributed pseudo-random numbers and
  42 //      functions that use those generators to generate numbers with
  43 //      different distributions:
  44 //
  45 //      class Rand32
  46 //      class Rand48
  47 //      solidSphereRand()
  48 //      hollowSphereRand()
  49 //      gaussRand()
  50 //      gaussSphereRand()
  51 //
  52 //-----------------------------------------------------------------------------
  53
  54 //
  55 // Here is the copyright for the *rand48() functions implemented for
  56 // Windows.
  57 //
  58
  59 //
  60 // Copyright (c) 1993 Martin Birgmeier
  61 // All rights reserved.
  62 //
  63 // You may redistribute unmodified or modified versions of this source
  64 // code provided that the above copyright notice and this and the
  65 // following conditions are retained.
  66 //
  67 // This software is provided ``as is'', and comes with no warranties
  68 // of any kind. I shall in no event be liable for anything that happens
  69 // to anyone/anything when using this software.
  70 //
  71
  72 #include <stdlib.h>
  73 #include <math.h>
  74
  75 namespace Imath {
  76
  77
  78 //-----------------------------------------------
  79 // Fast random-number generator that generates
  80 // a uniformly distributed sequence with a period
  81 // length of 2^32.
  82 //-----------------------------------------------
  83
  84 class Rand32
  85 {
  86   public:
  87
  88     //------------
  89     // Constructor
  90     //------------
  91
  92     Rand32 (unsigned long int seed = 0);
  93
  94
  95     //--------------------------------
  96     // Re-initialize with a given seed
  97     //--------------------------------
  98
  99     void                init (unsigned long int seed);
 100
 101
 102     //----------------------------------------------------------
 103     // Get the next value in the sequence (range: [false, true])
 104     //----------------------------------------------------------
 105
 106     bool                nextb ();
 107
 108
 109     //---------------------------------------------------------------
 110     // Get the next value in the sequence (range: [0 ... 0xffffffff])
 111     //---------------------------------------------------------------
 112
 113     unsigned long int   nexti ();
 114
 115
 116     //------------------------------------------------------
 117     // Get the next value in the sequence (range: [0 ... 1[)
 118     //------------------------------------------------------
 119
 120     float               nextf ();
 121
 122
 123     //-------------------------------------------------------------------
 124     // Get the next value in the sequence (range [rangeMin ... rangeMax[)
 125     //-------------------------------------------------------------------
 126
 127     float               nextf (float rangeMin, float rangeMax);
 128
 129
 130   private:
 131
 132     void                next ();
 133
 134     unsigned long int   _state;
 135 };
 136
 137
 138 //--------------------------------------------------------
 139 // Random-number generator based on the C Standard Library
 140 // functions drand48(), lrand48() & company; generates a
 141 // uniformly distributed sequence.
 142 //--------------------------------------------------------
 143
 144 class Rand48
 145 {
 146   public:
 147
 148     //------------
 149     // Constructor
 150     //------------
 151
 152     Rand48 (unsigned long int seed = 0);
 153
 154
 155     //--------------------------------
 156     // Re-initialize with a given seed
 157     //--------------------------------
 158
 159     void                init (unsigned long int seed);
 160
 161
 162     //----------------------------------------------------------
 163     // Get the next value in the sequence (range: [false, true])
 164     //----------------------------------------------------------
 165
 166     bool                nextb ();
 167
 168
 169     //---------------------------------------------------------------
 170     // Get the next value in the sequence (range: [0 ... 0x7fffffff])
 171     //---------------------------------------------------------------
 172
 173     long int            nexti ();
 174
 175
 176     //------------------------------------------------------
 177     // Get the next value in the sequence (range: [0 ... 1[)
 178     //------------------------------------------------------
 179
 180     double              nextf ();
 181
 182
 183     //-------------------------------------------------------------------
 184     // Get the next value in the sequence (range [rangeMin ... rangeMax[)
 185     //-------------------------------------------------------------------
 186
 187     double              nextf (double rangeMin, double rangeMax);
 188
 189
 190   private:
 191
 192     unsigned short int  _state[3];
 193
 194 #if defined ( _WIN32 ) || defined ( _WIN64 ) || defined ( __MWERKS__ )
 195     void                shiftState();
 196 #endif
 197 };
 198
 199
 200 //------------------------------------------------------------
 201 // Return random points uniformly distributed in a sphere with
 202 // radius 1 around the origin (distance from origin <= 1).
 203 //------------------------------------------------------------
 204
 205 template <class Vec, class Rand>
 206 Vec
 207 solidSphereRand (Rand &rand);
 208
 209
 210 //-------------------------------------------------------------
 211 // Return random points uniformly distributed on the surface of
 212 // a sphere with radius 1 around the origin.
 213 //-------------------------------------------------------------
 214
 215 template <class Vec, class Rand>
 216 Vec
 217 hollowSphereRand (Rand &rand);
 218
 219
 220 //-----------------------------------------------
 221 // Return random numbers with a normal (Gaussian)
 222 // distribution with zero mean and unit variance.
 223 //-----------------------------------------------
 224
 225 template <class Rand>
 226 float
 227 gaussRand (Rand &rand);
 228
 229
 230 //----------------------------------------------------
 231 // Return random points whose distance from the origin
 232 // has a normal (Gaussian) distribution with zero mean
 233 // and unit variance.
 234 //----------------------------------------------------
 235
 236 template <class Vec, class Rand>
 237 Vec
 238 gaussSphereRand (Rand &rand);
 239
 240
 241 //---------------
 242 // Implementation
 243 //---------------
 244
 245
 246 inline void
 247 Rand32::init (unsigned long int seed)
 248 {
 249     _state = (seed * 0xa5a573a5L) ^ 0x5a5a5a5aL;
 250 }
 251
 252
 253 inline
 254 Rand32::Rand32 (unsigned long int seed)
 255 {
 256     init (seed);
 257 }
 258
 259
 260 inline void
 261 Rand32::next ()
 262 {
 263     _state = 1664525L * _state + 1013904223L;
 264 }
 265
 266
 267 inline bool
 268 Rand32::nextb ()
 269 {
 270     next ();
 271     // Return the 31st (most significant) bit, by and-ing with 2 ^ 31.
 272     return !!(_state & 2147483648UL);
 273 }
 274
 275
 276 inline unsigned long int
 277 Rand32::nexti ()
 278 {
 279     next ();
 280     return _state & 0xffffffff;
 281 }
 282
 283
 284 inline float
 285 Rand32::nextf ()
 286 {
 287     next ();
 288     return ((int) (_state & 0xffffff)) * ((float) (1.0F / 0x1000000));
 289 }
 290
 291
 292 inline float
 293 Rand32::nextf (float rangeMin, float rangeMax)
 294 {
 295     return rangeMin + nextf() * (rangeMax - rangeMin);
 296 }
 297
 298
 299 inline void
 300 Rand48::init (unsigned long int seed)
 301 {
 302     seed = (seed * 0xa5a573a5L) ^ 0x5a5a5a5aL;
 303
 304     _state[0] = (unsigned short int) (seed);
 305     _state[1] = (unsigned short int) (seed >> 16);
 306     _state[2] = (unsigned short int) (seed);
 307 }
 308
 309
 310 inline
 311 Rand48::Rand48 (unsigned long int seed)
 312 {
 313     init (seed);
 314 }
 315
 316
 317 #if defined ( _WIN32 ) || defined ( _WIN64 ) || defined ( __MWERKS__ )
 318
 319 inline void
 320 Rand48::shiftState()
 321 {
 322     unsigned long   accu;
 323     unsigned short  temp[2];
 324
 325     accu = 0xe66dUL * ( unsigned long )_state[0] + 0x000bUL;
 326
 327     temp[0] = ( unsigned short )accu;   /* lower 16 bits */
 328     accu >>= sizeof( unsigned short ) * 8;
 329
 330     accu += 0xe66dUL * ( unsigned long )_state[1] +
 331             0xdeecUL * ( unsigned long )_state[0];
 332
 333     temp[1] = ( unsigned short )accu;   /* middle 16 bits */
 334     accu >>= sizeof( unsigned short ) * 8;
 335
 336     accu += 0xe66dUL * _state[2] +
 337             0xdeecUL * _state[1] +
 338             0x0005UL * _state[0];
 339
 340     _state[0] = temp[0];
 341     _state[1] = temp[1];
 342     _state[2] = ( unsigned short )accu;
 343 }
 344
 345 #endif
 346
 347 inline bool
 348 Rand48::nextb ()
 349 {
 350 #if defined ( _WIN32 ) || defined ( _WIN64 ) || defined ( __MWERKS__ )
 351     shiftState();
 352     return ( ( long( _state[2] ) << 15 ) + ( long( _state[1] ) >> 1 ) ) & 0x1;
 353 #else
 354     return nrand48 (_state) & 1;
 355 #endif
 356 }
 357
 358
 359 inline long int
 360 Rand48::nexti ()
 361 {
 362 #if defined ( _WIN32 ) || defined ( _WIN64 ) || defined ( __MWERKS__ )
 363     shiftState();
 364     return ( long( _state[2] ) << 15 ) + ( long( _state[1] ) >> 1 );
 365 #else
 366     return nrand48 (_state);
 367 #endif
 368 }
 369
 370
 371 inline double
 372 Rand48::nextf ()
 373 {
 374 #if defined ( _WIN32 ) || defined ( _WIN64 ) || defined ( __MWERKS__ )
 375     shiftState();
 376     return ldexp( double( _state[0] ), -48 ) +
 377            ldexp( double( _state[1] ), -32 ) +
 378            ldexp( double( _state[2] ), -16 );
 379 #else
 380     return erand48 (_state);
 381 #endif
 382 }
 383
 384
 385 inline double
 386 Rand48::nextf (double rangeMin, double rangeMax)
 387 {
 388     return rangeMin + nextf() * (rangeMax - rangeMin);
 389 }
 390
 391
 392 template <class Vec, class Rand>
 393 Vec
 394 solidSphereRand (Rand &rand)
 395 {
 396     Vec v;
 397
 398     do
 399     {
 400         for (unsigned int i = 0; i < Vec::dimensions(); i++)
 401             v[i] = (typename Vec::BaseType) rand.nextf (-1, 1);
 402     }
 403     while (v.length2() > 1);
 404
 405     return v;
 406 }
 407
 408
 409 template <class Vec, class Rand>
 410 Vec
 411 hollowSphereRand (Rand &rand)
 412 {
 413     Vec v;
 414     typename Vec::BaseType length;
 415
 416     do
 417     {
 418         for (unsigned int i = 0; i < Vec::dimensions(); i++)
 419             v[i] = (typename Vec::BaseType) rand.nextf (-1, 1);
 420
 421         length = v.length();
 422     }
 423     while (length > 1 || length == 0);
 424
 425     return v / length;
 426 }
 427
 428
 429 template <class Rand>
 430 float
 431 gaussRand (Rand &rand)
 432 {
 433     float x;            // Note: to avoid numerical problems with very small
 434     float y;            // numbers, we make these variables singe-precision
 435     float length2;      // floats, but later we call the double-precision log()
 436                         // and sqrt() functions instead of logf() and sqrtf().
 437     do
 438     {
 439         x = float (rand.nextf (-1, 1));
 440         y = float (rand.nextf (-1, 1));
 441         length2 = x * x + y * y;
 442     }
 443     while (length2 >= 1 || length2 == 0);
 444
 445     return x * sqrt (-2 * log (length2) / length2);
 446 }
 447
 448
 449 template <class Vec, class Rand>
 450 Vec
 451 gaussSphereRand (Rand &rand)
 452 {
 453     return hollowSphereRand <Vec> (rand) * gaussRand (rand);
 454 }
 455
 456 double drand48();
 457 long int lrand48();
 458
 459 } // namespace Imath
 460
 461 #endif