1 /*******************************************************************************
2 * randomsequences.cpp
3 *
4 * ---------------------------------------------------------------------------
5 * Persistence of Vision Ray Tracer ('POV-Ray') version 3.7.
6 * Copyright 1991-2013 Persistence of Vision Raytracer Pty. Ltd.
7 *
8 * POV-Ray is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU Affero General Public License as
10 * published by the Free Software Foundation, either version 3 of the
11 * License, or (at your option) any later version.
12 *
13 * POV-Ray is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU Affero General Public License for more details.
17 *
18 * You should have received a copy of the GNU Affero General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 * ---------------------------------------------------------------------------
21 * POV-Ray is based on the popular DKB raytracer version 2.12.
22 * DKBTrace was originally written by David K. Buck.
23 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
24 * ---------------------------------------------------------------------------
25 * $File: //depot/povray/smp/source/backend/support/randomsequences.cpp $
26 * $Revision: #23 $
27 * $Change: 6132 $
28 * $DateTime: 2013/11/25 14:23:41 $
29 * $Author: clipka $
30 *******************************************************************************/
31
32 #include <cassert>
33 #include <stdexcept>
34 #include <map>
35
36 #include <boost/random/mersenne_twister.hpp>
37 #include <boost/random/uniform_int.hpp>
38 #include <boost/random/uniform_real.hpp>
39 #include <boost/random/variate_generator.hpp>
40
41 // frame.h must always be the first POV file included (pulls in platform config)
42 #include "backend/frame.h"
43 #include "backend/support/randomsequences.h"
44
45 // this must be the last file included
46 #include "base/povdebug.h"
47
48 namespace pov
49 {
50
51 using namespace pov_base;
52
53 using boost::uniform_int;
54 using boost::uniform_real;
55 using boost::variate_generator;
56 using boost::mt19937;
57
58 #ifndef SIZE_MAX
59 #define SIZE_MAX ((size_t)-1)
60 #endif
61
62 #define PRIME_TABLE_COUNT 25
63 unsigned int primeTable[PRIME_TABLE_COUNT] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
64
65
66
67 /*****************************************************************************
68 *
69 * FUNCTION
70 *
71 * stream_rand
72 *
73 * INPUT
74 *
75 * stream - number of random stream
76 *
77 * OUTPUT
78 *
79 * RETURNS
80 *
81 * DBL - random value
82 *
83 * AUTHOR
84 *
85 * Dieter Bayer
86 *
87 * DESCRIPTION
88 *
89 * Standard pseudo-random function.
90 *
91 * CHANGES
92 *
93 * Feb 1996 : Creation.
94 * Mar 1996 : Return 2^32 random values instead of 2^16 [AED]
95 *
96 ******************************************************************************/
97
POV_rand(unsigned int & next_rand)98 DBL POV_rand(unsigned int& next_rand)
99 {
100 next_rand = next_rand * 1812433253L + 12345L;
101
102 return((DBL)(next_rand & 0xFFFFFFFFUL) / 0xFFFFFFFFUL);
103 }
104
105
106 /**********************************************************************************
107 * Legacy Code
108 *********************************************************************************/
109
RandomInts(int minval,int maxval,size_t count)110 vector<int> RandomInts(int minval, int maxval, size_t count)
111 {
112 mt19937 generator;
113 uniform_int<int> distribution(minval, maxval);
114 variate_generator<mt19937, uniform_int<int> > sequence(generator, distribution);
115 vector<int> rands(count);
116
117 for(size_t i = 0; i < count; i++)
118 rands[i] = sequence();
119
120 return rands;
121 }
122
RandomDoubles(double minval,double maxval,size_t count)123 vector<double> RandomDoubles(double minval, double maxval, size_t count)
124 {
125 mt19937 generator;
126 uniform_real<double> distribution(minval, maxval);
127 variate_generator<mt19937, uniform_real<double> > sequence(generator, distribution);
128 vector<double> rands(count);
129
130 for(size_t i = 0; i < count; i++)
131 rands[i] = sequence();
132
133 return rands;
134 }
135
RandomIntSequence(int minval,int maxval,size_t count)136 RandomIntSequence::RandomIntSequence(int minval, int maxval, size_t count) :
137 values(RandomInts(minval, maxval, count))
138 {
139 }
140
Generator(RandomIntSequence * seq,size_t seedindex)141 RandomIntSequence::Generator::Generator(RandomIntSequence *seq, size_t seedindex) :
142 sequence(seq),
143 index(seedindex)
144 {
145 }
146
operator ()(size_t seedindex)147 int RandomIntSequence::operator()(size_t seedindex)
148 {
149 seedindex = seedindex % values.size();
150 return values[seedindex];
151 }
152
operator ()()153 int RandomIntSequence::Generator::operator()()
154 {
155 index = (index + 1) % sequence->values.size();
156 return (*sequence)(index);
157 }
158
operator ()(size_t seedindex)159 int RandomIntSequence::Generator::operator()(size_t seedindex)
160 {
161 return (*sequence)(seedindex);
162 }
163
GetSeed() const164 size_t RandomIntSequence::Generator::GetSeed() const
165 {
166 return index;
167 }
168
SetSeed(size_t seedindex)169 void RandomIntSequence::Generator::SetSeed(size_t seedindex)
170 {
171 index = seedindex % sequence->values.size();
172 }
173
RandomDoubleSequence(double minval,double maxval,size_t count)174 RandomDoubleSequence::RandomDoubleSequence(double minval, double maxval, size_t count) :
175 values(RandomDoubles(minval, maxval, count))
176 {
177 }
178
Generator(RandomDoubleSequence * seq,size_t seedindex)179 RandomDoubleSequence::Generator::Generator(RandomDoubleSequence *seq, size_t seedindex) :
180 sequence(seq),
181 index(seedindex)
182 {
183 }
184
operator ()(size_t seedindex)185 double RandomDoubleSequence::operator()(size_t seedindex)
186 {
187 seedindex = seedindex % values.size();
188 return values[seedindex];
189 }
190
operator ()()191 double RandomDoubleSequence::Generator::operator()()
192 {
193 index = (index + 1) % sequence->values.size();
194 return (*sequence)(index);
195 }
196
operator ()(size_t seedindex)197 double RandomDoubleSequence::Generator::operator()(size_t seedindex)
198 {
199 return (*sequence)(seedindex);
200 }
201
GetSeed() const202 size_t RandomDoubleSequence::Generator::GetSeed() const
203 {
204 return index;
205 }
206
SetSeed(size_t seedindex)207 void RandomDoubleSequence::Generator::SetSeed(size_t seedindex)
208 {
209 index = seedindex % sequence->values.size();
210 }
211
212
213 /**********************************************************************************
214 * Local Types : Abstract Generators
215 *********************************************************************************/
216
217 /**
218 * Abstract template class representing a generator for numbers that can be accessed both sequentially and by index.
219 */
220 template<class Type>
221 class HybridNumberGenerator : public SeedableNumberGenerator<Type>, public IndexedNumberGenerator<Type>
222 {
223 public:
224
225 HybridNumberGenerator(size_t size = 0);
226 virtual Type operator()();
227 virtual shared_ptr<vector<Type> > GetSequence(size_t count);
228 virtual size_t MaxIndex() const;
229 virtual size_t CycleLength() const;
230 virtual void Seed(size_t seed);
231
232 protected:
233
234 const size_t size;
235 size_t index;
236 };
237
238
239 /**********************************************************************************
240 * Local Types : Linear Generators
241 *********************************************************************************/
242
243 /**
244 * Template class representing a generator for uniformly distributed numbers in a given range.
245 */
246 template<class Type, class BoostGenerator, class UniformType, size_t CYCLE_LENGTH = SIZE_MAX>
247 class UniformRandomNumberGenerator : public SequentialNumberGenerator<Type>
248 {
249 public:
250
251 struct ParameterStruct {
252 ParameterStruct(Type minval, Type maxval);
253 Type minval, maxval;
254 bool operator< (const ParameterStruct& other) const;
255 };
256
257 UniformRandomNumberGenerator(const ParameterStruct& param);
258 UniformRandomNumberGenerator(Type minval, Type maxval);
259 virtual Type operator()();
260 virtual size_t CycleLength() const;
261
262 protected:
263 variate_generator<BoostGenerator, UniformType> generator;
264 };
265
266 typedef UniformRandomNumberGenerator<int, mt19937, uniform_int<int> > Mt19937IntGenerator;
267 typedef UniformRandomNumberGenerator<double, mt19937, uniform_real<double> > Mt19937DoubleGenerator;
268
269 /**
270 * Generator for a 1-dimensional Halton sequence (aka van-der-Corput sequence).
271 * This class fulfills the boost UniformRandomNumberGenerator requirements,
272 * except that the numbers generated are actually sub-random.
273 */
274 template<class Type>
275 class HaltonGenerator : public HybridNumberGenerator<Type>
276 {
277 public:
278
279 struct ParameterStruct {
280 ParameterStruct(unsigned int base, Type minval, Type maxval);
281 unsigned int base;
282 Type minval, maxval;
283 bool operator< (const ParameterStruct& other) const;
284 };
285
286 HaltonGenerator(const ParameterStruct& param);
287 HaltonGenerator(unsigned int base, Type minval, Type maxval);
288 /// Returns a particular number from the sequence.
289 virtual double operator[](size_t index) const;
290
291 protected:
292
293 unsigned int base;
294 Type minval;
295 Type scale;
296 };
297
298 typedef HaltonGenerator<int> HaltonIntGenerator;
299 typedef HaltonGenerator<double> HaltonDoubleGenerator;
300
301
302 /**********************************************************************************
303 * Local Types : Vector Generators
304 *********************************************************************************/
305
306 /**
307 * Class generating a cosine-weighted hemispherical direction vector compatible with earlier POV-Ray versions.
308 * This class uses a 1600-element hard-coded directions originally used for radiosity.
309 */
310 class LegacyCosWeightedDirectionGenerator : public HybridNumberGenerator<Vector3d>
311 {
312 public:
313
314 static const int NumEntries = 1600;
315
316 struct ParameterStruct
317 {
318 bool operator< (const ParameterStruct& other) const;
319 };
320
321 LegacyCosWeightedDirectionGenerator(const ParameterStruct& dummy);
322 virtual Vector3d operator[](size_t i) const;
323 };
324
325
326 /**
327 * Abstract template class generating a vector based on a 2D Halton sequence.
328 */
329 template<class Type, class TypeA, class TypeB = TypeA>
330 class Halton2dBasedGenerator : public HybridNumberGenerator<Type>
331 {
332 public:
333
334 struct ParameterStruct
335 {
336 ParameterStruct(unsigned int baseA, unsigned int baseB, TypeA minvalA, TypeA maxvalA, TypeB minvalB, TypeB maxvalB);
337 unsigned int baseA, baseB;
338 TypeA minvalA, maxvalA;
339 TypeB minvalB, maxvalB;
340 bool operator< (const ParameterStruct& other) const;
341 };
342
343 Halton2dBasedGenerator(const ParameterStruct& param);
344 virtual Type operator[](size_t i) const = 0;
345
346 protected:
347
348 shared_ptr<HaltonDoubleGenerator> generatorA;
349 shared_ptr<HaltonDoubleGenerator> generatorB;
350 };
351
352 /**
353 * Class generating cosine-weighted hemispherical direction vectors, centered around the Y axis, based on a 2D Halton sequence.
354 */
355 class HaltonCosWeightedDirectionGenerator : public Halton2dBasedGenerator<Vector3d, double>
356 {
357 public:
358
359 struct ParameterStruct : public Halton2dBasedGenerator<Vector3d, double>::ParameterStruct
360 {
361 ParameterStruct(unsigned int baseA, unsigned int baseB);
362 };
363
364 HaltonCosWeightedDirectionGenerator(const ParameterStruct& param);
365 virtual Vector3d operator[](size_t i) const;
366 };
367
368 /**
369 * Class generating uniformly distributed points within the unit circle based on a 2D Halton sequence.
370 */
371 class HaltonOnDiscGenerator : public Halton2dBasedGenerator<Vector2d, double>
372 {
373 public:
374
375 struct ParameterStruct : public Halton2dBasedGenerator<Vector2d, double>::ParameterStruct
376 {
377 ParameterStruct(unsigned int baseA, unsigned int baseB, double radius);
378 };
379
380 HaltonOnDiscGenerator(const ParameterStruct& param);
381 virtual Vector2d operator[](size_t i) const;
382 };
383
384 /**
385 * Class generating uniformly distributed points on the unit sphere based on a 2D Halton sequence.
386 */
387 class HaltonUniformDirectionGenerator : public Halton2dBasedGenerator<Vector3d, double>
388 {
389 public:
390
391 struct ParameterStruct : public Halton2dBasedGenerator<Vector3d, double>::ParameterStruct
392 {
393 ParameterStruct(unsigned int baseA, unsigned int baseB);
394 };
395
396 HaltonUniformDirectionGenerator(const ParameterStruct& param);
397 virtual Vector3d operator[](size_t i) const;
398 };
399
400 /**
401 * Class generating uniformly distributed points within a square based on a 2D Halton sequence.
402 */
403 class Halton2dGenerator : public Halton2dBasedGenerator<Vector2d, double>
404 {
405 public:
406 Halton2dGenerator(const ParameterStruct& param);
407 virtual Vector2d operator[](size_t i) const;
408 };
409
410
411 /**********************************************************************************
412 * Local Types : Auxiliary
413 *********************************************************************************/
414
415 /**
416 * Template class representing a factory for pre-computed number tables.
417 */
418 template<class Type>
419 class NumberSequenceFactory
420 {
421 public:
422
423 /// Sets up the factory to use a given sequence.
424 NumberSequenceFactory(shared_ptr<vector<Type> const> masterSequence);
425 /// Sets up the factory to use a given number source.
426 NumberSequenceFactory(shared_ptr<SequentialNumberGenerator<Type> > master);
427 /// Sets up the factory to use a given number source, pre-computing a given number of elements.
428 NumberSequenceFactory(shared_ptr<SequentialNumberGenerator<Type> > master, size_t count);
429 /// Gets a reference to a table of pre-computed numbers having at least the given size.
430 /// @note The vector returned may contain more elements than requested.
431 shared_ptr<vector<Type> const> operator()(size_t count);
432
433 protected:
434
435 typedef SequentialNumberGenerator<Type> Generator;
436 typedef shared_ptr<Generator> GeneratorPtr;
437 typedef vector<Type> Sequence;
438 typedef shared_ptr<Sequence> SequencePtr;
439 typedef shared_ptr<Sequence const> SequenceConstPtr;
440
441 GeneratorPtr master;
442 SequenceConstPtr masterSequence;
443 boost::mutex masterMutex;
444 };
445
446 typedef NumberSequenceFactory<int> IntSequenceFactory;
447 typedef NumberSequenceFactory<double> DoubleSequenceFactory;
448 typedef NumberSequenceFactory<Vector3d> VectorSequenceFactory;
449
450
451 /**
452 * Template class representing a meta-factory for factories for pre-computed number tables.
453 */
454 template<class ValueType, class GeneratorType>
455 class NumberSequenceMetaFactory
456 {
457 public:
458
459 static shared_ptr<NumberSequenceFactory<ValueType> > GetFactory(const typename GeneratorType::ParameterStruct& param);
460
461 protected:
462
463 typedef NumberSequenceFactory<ValueType> Factory;
464 typedef shared_ptr<Factory> FactoryPtr;
465 typedef weak_ptr<Factory> FactoryWeakPtr;
466 typedef std::map<typename GeneratorType::ParameterStruct, FactoryWeakPtr> FactoryTable;
467
468 static FactoryTable* lookupTable;
469 static boost::mutex lookupMutex;
470 };
471
472 typedef NumberSequenceMetaFactory<int, Mt19937IntGenerator> Mt19937IntMetaFactory;
473 typedef NumberSequenceMetaFactory<double, Mt19937DoubleGenerator> Mt19937DoubleMetaFactory;
474 typedef NumberSequenceMetaFactory<Vector3d, LegacyCosWeightedDirectionGenerator> LegacyCosWeightedDirectionMetaFactory;
475 typedef NumberSequenceMetaFactory<Vector3d, HaltonCosWeightedDirectionGenerator> HaltonCosWeightedDirectionMetaFactory;
476 typedef NumberSequenceMetaFactory<double, HaltonDoubleGenerator> HaltonUniformDoubleMetaFactory;
477 typedef NumberSequenceMetaFactory<Vector3d, HaltonUniformDirectionGenerator> HaltonUniformDirectionMetaFactory;
478 typedef NumberSequenceMetaFactory<Vector2d, HaltonOnDiscGenerator> HaltonOnDiscMetaFactory;
479 typedef NumberSequenceMetaFactory<Vector2d, Halton2dGenerator> Halton2dMetaFactory;
480
481
482 /**
483 * Template class representing a generator for pre-computed numbers using a shared values table.
484 */
485 template<class Type>
486 class PrecomputedNumberGenerator : public HybridNumberGenerator<Type>
487 {
488 public:
489
490 /// Construct from a sequence factory.
PrecomputedNumberGenerator(shared_ptr<NumberSequenceFactory<Type>> master,size_t size)491 PrecomputedNumberGenerator(shared_ptr<NumberSequenceFactory<Type> > master, size_t size) :
492 HybridNumberGenerator<Type>(size),
493 values((*master)(size))
494 {}
495
496 /// Returns a particular number from the sequence.
operator [](size_t i) const497 virtual Type operator[](size_t i) const
498 {
499 // According to C++ standard, template classes cannot refer to parent template classes' members by unqualified name
500 const size_t& size = HybridNumberGenerator<Type>::size;
501 return (*values)[i % size];
502 }
503 /// Returns a particular subset from the sequence.
GetSequence(size_t index,size_t count) const504 virtual shared_ptr<vector<Type> > GetSequence(size_t index, size_t count) const
505 {
506 // According to C++ standard, template classes cannot refer to parent template classes' members by unqualified name
507 const size_t& size = HybridNumberGenerator<Type>::size;
508 shared_ptr<vector<Type> > data(new vector<Type>);
509 data->reserve(count);
510 size_t i = index % size;
511 while (count >= size - i) // handle wrap-around
512 {
513 data->insert(data->end(), values->begin() + i, values->begin() + size);
514 count -= (size - i);
515 i = 0;
516 }
517 data->insert(data->end(), values->begin() + i, values->begin() + i + count);
518 return data;
519 }
520
521 protected:
522
523 shared_ptr<vector<Type> const> values;
524 };
525
526 typedef PrecomputedNumberGenerator<int> PrecomputedIntGenerator;
527 typedef PrecomputedNumberGenerator<double> PrecomputedDoubleGenerator;
528 typedef PrecomputedNumberGenerator<Vector3d> PrecomputedVectorGenerator;
529 typedef PrecomputedNumberGenerator<Vector2d> PrecomputedVector2dGenerator;
530
531
532 /**********************************************************************************
533 * HybridNumberGenerator implementation
534 *********************************************************************************/
535
536 template<class Type>
HybridNumberGenerator(size_t size)537 HybridNumberGenerator<Type>::HybridNumberGenerator(size_t size) :
538 size(size),
539 index(0)
540 {}
541
542 template<class Type>
operator ()()543 Type HybridNumberGenerator<Type>::operator()()
544 {
545 const Type& data = (*this)[index ++];
546 if (size != 0)
547 index = index % size;
548 return data;
549 }
550
551 template<class Type>
GetSequence(size_t count)552 shared_ptr<vector<Type> > HybridNumberGenerator<Type>::GetSequence(size_t count)
553 {
554 shared_ptr<vector<Type> > data(IndexedNumberGenerator<Type>::GetSequence(index, count));
555 index += count;
556 if (size != 0)
557 index = index % size;
558 return data;
559 }
560
561 template<class Type>
MaxIndex() const562 size_t HybridNumberGenerator<Type>::MaxIndex() const
563 {
564 return size - 1;
565 }
566
567 template<class Type>
CycleLength() const568 size_t HybridNumberGenerator<Type>::CycleLength() const
569 {
570 return size;
571 }
572
573 template<class Type>
Seed(size_t seed)574 void HybridNumberGenerator<Type>::Seed(size_t seed)
575 {
576 index = seed % size;
577 }
578
579
580 /**********************************************************************************
581 * UniformRandomNumberGenerator implementation
582 *********************************************************************************/
583
584 template<class Type, class BoostGenerator, class UniformType, size_t CYCLE_LENGTH>
ParameterStruct(Type minval,Type maxval)585 UniformRandomNumberGenerator<Type,BoostGenerator,UniformType,CYCLE_LENGTH>::ParameterStruct::ParameterStruct(Type minval, Type maxval) :
586 minval(minval), maxval(maxval)
587 {}
588
589 template<class Type, class BoostGenerator, class UniformType, size_t CYCLE_LENGTH>
operator <(const ParameterStruct & other) const590 bool UniformRandomNumberGenerator<Type,BoostGenerator,UniformType,CYCLE_LENGTH>::ParameterStruct::operator< (const ParameterStruct& other) const
591 {
592 if (minval != other.minval)
593 return (minval < other.minval);
594 else
595 return (maxval < other.maxval);
596 }
597
598 template<class Type, class BoostGenerator, class UniformType, size_t CYCLE_LENGTH>
UniformRandomNumberGenerator(const ParameterStruct & param)599 UniformRandomNumberGenerator<Type,BoostGenerator,UniformType,CYCLE_LENGTH>::UniformRandomNumberGenerator(const ParameterStruct& param) :
600 generator(BoostGenerator(), UniformType(param.minval, param.maxval))
601 {}
602
603 template<class Type, class BoostGenerator, class UniformType, size_t CYCLE_LENGTH>
UniformRandomNumberGenerator(Type minval,Type maxval)604 UniformRandomNumberGenerator<Type,BoostGenerator,UniformType,CYCLE_LENGTH>::UniformRandomNumberGenerator(Type minval, Type maxval) :
605 generator(BoostGenerator(), UniformType(minval, maxval))
606 {}
607
608 template<class Type, class BoostGenerator, class UniformType, size_t CYCLE_LENGTH>
609 Type UniformRandomNumberGenerator<Type,BoostGenerator,UniformType,CYCLE_LENGTH>::operator()()
610 {
611 return generator();
612 }
613
614 template<class Type, class BoostGenerator, class UniformType, size_t CYCLE_LENGTH>
CycleLength() const615 size_t UniformRandomNumberGenerator<Type,BoostGenerator,UniformType,CYCLE_LENGTH>::CycleLength() const
616 {
617 return CYCLE_LENGTH;
618 }
619
620
621 /**********************************************************************************
622 * HaltonGenerator implementation
623 *********************************************************************************/
624
625 template<class Type>
ParameterStruct(unsigned int base,Type minval,Type maxval)626 HaltonGenerator<Type>::ParameterStruct::ParameterStruct(unsigned int base, Type minval, Type maxval) :
627 base(base), minval(minval), maxval(maxval)
628 {}
629
630 template<class Type>
operator <(const ParameterStruct & other) const631 bool HaltonGenerator<Type>::ParameterStruct::operator< (const ParameterStruct& other) const
632 {
633 if (base != other.base)
634 return (base < other.base);
635 else if (minval != other.minval)
636 return (minval < other.minval);
637 else
638 return (maxval < other.maxval);
639 }
640
641 template<class Type>
HaltonGenerator(const ParameterStruct & param)642 HaltonGenerator<Type>::HaltonGenerator(const ParameterStruct& param) :
643 base(param.base),
644 minval(param.minval),
645 scale(param.maxval-param.minval)
646 {
647 }
648
649 template<class Type>
HaltonGenerator(unsigned int base,Type minval,Type maxval)650 HaltonGenerator<Type>::HaltonGenerator(unsigned int base, Type minval, Type maxval) :
651 base(base),
652 minval(minval),
653 scale(maxval-minval)
654 {
655 }
656
657 template<class Type>
operator [](size_t index) const658 double HaltonGenerator<Type>::operator[](size_t index) const
659 {
660 size_t i = 1 + index; // index starts at 0, but halton sequence as implemented here starts at 1
661
662 double h = 0;
663 double q = 1.0/base;
664 unsigned int digit;
665
666 while (i > 0)
667 {
668 digit = (unsigned int)(i % base);
669 h = h + digit * q;
670 i /= base;
671 q /= base;
672 }
673
674 return minval + (Type)(h * scale);
675 }
676
677
678 /**********************************************************************************
679 * NumberSequenceFactory implementation
680 *********************************************************************************/
681
682 template<class Type>
NumberSequenceFactory(shared_ptr<vector<Type> const> masterSequence)683 NumberSequenceFactory<Type>::NumberSequenceFactory(shared_ptr<vector<Type> const> masterSequence) :
684 masterSequence(masterSequence)
685 {}
686
687 template<class Type>
NumberSequenceFactory(shared_ptr<SequentialNumberGenerator<Type>> master)688 NumberSequenceFactory<Type>::NumberSequenceFactory(shared_ptr<SequentialNumberGenerator<Type> > master) :
689 master(master)
690 {}
691
692 template<class Type>
NumberSequenceFactory(shared_ptr<SequentialNumberGenerator<Type>> master,size_t count)693 NumberSequenceFactory<Type>::NumberSequenceFactory(shared_ptr<SequentialNumberGenerator<Type> > master, size_t count) :
694 master(master)
695 {
696 (*this)(count); // force initial sequence to be generated
697 }
698
699 template<class Type>
operator ()(size_t count)700 shared_ptr<vector<Type> const> NumberSequenceFactory<Type>::operator()(size_t count)
701 {
702 boost::mutex::scoped_lock lock(masterMutex);
703 if (!masterSequence)
704 {
705 // No values pre-computed yet; do it now.
706 masterSequence = SequenceConstPtr(master->GetSequence(count));
707 }
708 else if ((masterSequence->size() < count) && master)
709 {
710 // Not enough values pre-computed; release the current values list and build a larger one.
711 // NB: We're not simply appending to the current values list, because that might require re-allocation
712 // and interfere with other threads trying to read from the list. To avoid having to synchronize
713 // all read accesses, we're going for the less memory-efficient approach.
714 size_t newCount = count;
715 if (masterSequence->size() > newCount / 2)
716 {
717 // make sure to pre-compute at least twice the already-computed size, so we don't waste too much space with
718 if (masterSequence->size() > SIZE_MAX / 2) // play it safe (though that'll have us run out of memory anyway)
719 newCount = SIZE_MAX;
720 else
721 newCount = masterSequence->size() * 2;
722 }
723 // Pull more data from our master generator.
724 // NB: We're using a temporary pointer to the new values list, so we can keep the master list const,
725 // lest anyone might accidently modify it while other threads are reading it.
726 SequenceConstPtr newSequence(master->GetSequence(newCount - masterSequence->size()));
727 SequencePtr mergedSequence(new Sequence(*masterSequence));
728 mergedSequence->insert(mergedSequence->end(), newSequence->begin(), newSequence->end());
729 masterSequence = mergedSequence;
730 }
731 return masterSequence;
732 }
733
734
735 /**********************************************************************************
736 * NumberSequenceMetaFactory implementation
737 *********************************************************************************/
738
739 template<class ValueType, class GeneratorType>
740 std::map<typename GeneratorType::ParameterStruct, weak_ptr<NumberSequenceFactory<ValueType> > >* NumberSequenceMetaFactory<ValueType, GeneratorType>::lookupTable;
741
742 template<class ValueType, class GeneratorType>
743 boost::mutex NumberSequenceMetaFactory<ValueType, GeneratorType>::lookupMutex;
744
745 template<class ValueType, class GeneratorType>
GetFactory(const typename GeneratorType::ParameterStruct & param)746 shared_ptr<NumberSequenceFactory<ValueType> > NumberSequenceMetaFactory<ValueType, GeneratorType>::GetFactory(const typename GeneratorType::ParameterStruct& param)
747 {
748 boost::mutex::scoped_lock lock(lookupMutex);
749 if (!lookupTable)
750 lookupTable = new FactoryTable();
751 FactoryPtr factory = (*lookupTable)[param].lock();
752 if (!factory)
753 {
754 shared_ptr<GeneratorType> masterGenerator(new GeneratorType(param));
755 factory = FactoryPtr(new Factory(shared_ptr<SequentialNumberGenerator<ValueType> >(masterGenerator)));
756 (*lookupTable)[param] = factory;
757 }
758 return factory;
759 }
760
761
762 /**********************************************************************************
763 * LegacyCosWeightedDirectionGenerator implementation
764 *********************************************************************************/
765
766 extern BYTE_XYZ rad_samples[]; // defined in rad_data.cpp
767
operator <(const ParameterStruct & other) const768 bool LegacyCosWeightedDirectionGenerator::ParameterStruct::operator< (const ParameterStruct& other) const
769 {
770 return false; // all instances are equal
771 }
772
LegacyCosWeightedDirectionGenerator(const ParameterStruct & dummy)773 LegacyCosWeightedDirectionGenerator::LegacyCosWeightedDirectionGenerator(const ParameterStruct& dummy)
774 {}
775
operator [](size_t i) const776 Vector3d LegacyCosWeightedDirectionGenerator::operator[](size_t i) const
777 {
778 Vector3d result;
779 VUnpack(result, &(rad_samples[i % NumEntries]));
780 return result;
781 }
782
783
784 /**********************************************************************************
785 * Halton2dBasedGenerator implementation
786 *********************************************************************************/
787
788 template<class Type, class TypeA, class TypeB>
ParameterStruct(unsigned int baseA,unsigned int baseB,TypeA minvalA,TypeA maxvalA,TypeB minvalB,TypeB maxvalB)789 Halton2dBasedGenerator<Type, TypeA, TypeB>::ParameterStruct::ParameterStruct(unsigned int baseA, unsigned int baseB, TypeA minvalA, TypeA maxvalA, TypeB minvalB, TypeB maxvalB) :
790 baseA(baseA), baseB(baseB),
791 minvalA(minvalA), maxvalA(maxvalA),
792 minvalB(minvalB), maxvalB(maxvalB)
793 {}
794
795 template<class Type, class TypeA, class TypeB>
operator <(const ParameterStruct & other) const796 bool Halton2dBasedGenerator<Type, TypeA, TypeB>::ParameterStruct::operator< (const ParameterStruct& other) const
797 {
798 if (baseA != other.baseA)
799 return (baseA < other.baseA);
800 else if (baseB != other.baseB)
801 return (baseB < other.baseB);
802 else if (minvalA != other.minvalA)
803 return (minvalA < other.minvalA);
804 else if (maxvalA != other.maxvalA)
805 return (maxvalA < other.maxvalA);
806 else if (minvalB != other.minvalB)
807 return (minvalB < other.minvalB);
808 else
809 return (maxvalB < other.maxvalB);
810 }
811
812 template<class Type, class TypeA, class TypeB>
Halton2dBasedGenerator(const ParameterStruct & param)813 Halton2dBasedGenerator<Type, TypeA, TypeB>::Halton2dBasedGenerator(const ParameterStruct& param) :
814 generatorA(new HaltonDoubleGenerator(param.baseA, param.minvalA, param.maxvalA)),
815 generatorB(new HaltonDoubleGenerator(param.baseB, param.minvalB, param.maxvalB))
816 {}
817
818
819 /**********************************************************************************
820 * HaltonCosWeightedDirectionGenerator implementation
821 *********************************************************************************/
822
HaltonCosWeightedDirectionGenerator(const ParameterStruct & param)823 HaltonCosWeightedDirectionGenerator::HaltonCosWeightedDirectionGenerator(const ParameterStruct& param) :
824 Halton2dBasedGenerator<Vector3d,double,double>(param)
825 {}
826
ParameterStruct(unsigned int baseA,unsigned int baseB)827 HaltonCosWeightedDirectionGenerator::ParameterStruct::ParameterStruct(unsigned int baseA, unsigned int baseB) :
828 Halton2dBasedGenerator<Vector3d,double,double>::ParameterStruct(baseA, baseB, 0.0, 1.0, 0.0, 2*M_PI)
829 {}
830
operator [](size_t i) const831 Vector3d HaltonCosWeightedDirectionGenerator::operator[](size_t i) const
832 {
833 double r = sqrt((*generatorA)[i]);
834 double theta = (*generatorB)[i];
835 double x = r * cos(theta);
836 double z = r * sin(theta);
837 double y = sqrt (1 - x*x - z*z);
838 return Vector3d(x, y, z);
839 }
840
841
842 /**********************************************************************************
843 * HaltonOnDiscGenerator implementation
844 *********************************************************************************/
845
HaltonOnDiscGenerator(const ParameterStruct & param)846 HaltonOnDiscGenerator::HaltonOnDiscGenerator(const ParameterStruct& param) :
847 Halton2dBasedGenerator<Vector2d,double,double>(param)
848 {}
849
ParameterStruct(unsigned int baseA,unsigned int baseB,double radius)850 HaltonOnDiscGenerator::ParameterStruct::ParameterStruct(unsigned int baseA, unsigned int baseB, double radius) :
851 Halton2dBasedGenerator<Vector2d,double,double>::ParameterStruct(baseA, baseB, 0.0, radius*radius, 0.0, 2*M_PI)
852 {}
853
operator [](size_t i) const854 Vector2d HaltonOnDiscGenerator::operator[](size_t i) const
855 {
856 double r = sqrt((*generatorA)[i]);
857 double theta = (*generatorB)[i];
858 double x = r * cos(theta);
859 double y = r * sin(theta);
860 return Vector2d(x, y);
861 }
862
863
864 /**********************************************************************************
865 * Halton2dGenerator implementation
866 *********************************************************************************/
867
Halton2dGenerator(const ParameterStruct & param)868 Halton2dGenerator::Halton2dGenerator(const ParameterStruct& param) :
869 Halton2dBasedGenerator<Vector2d,double,double>(param)
870 {}
871
operator [](size_t i) const872 Vector2d Halton2dGenerator::operator[](size_t i) const
873 {
874 double x = (*generatorA)[i];
875 double y = (*generatorB)[i];
876 return Vector2d(x, y);
877 }
878
879
880 /**********************************************************************************
881 * HaltonUniformDirectionGenerator implementation
882 *********************************************************************************/
883
HaltonUniformDirectionGenerator(const ParameterStruct & param)884 HaltonUniformDirectionGenerator::HaltonUniformDirectionGenerator(const ParameterStruct& param) :
885 Halton2dBasedGenerator<Vector3d,double,double>(param)
886 {}
887
ParameterStruct(unsigned int baseA,unsigned int baseB)888 HaltonUniformDirectionGenerator::ParameterStruct::ParameterStruct(unsigned int baseA, unsigned int baseB) :
889 Halton2dBasedGenerator<Vector3d,double,double>::ParameterStruct(baseA, baseB, -1.0, 1.0, 0.0, 2*M_PI)
890 {}
891
operator [](size_t i) const892 Vector3d HaltonUniformDirectionGenerator::operator[](size_t i) const
893 {
894 double x = (*generatorA)[i];
895 double r = sqrt(1 - x*x);
896 double theta = (*generatorB)[i];
897 double y = r * cos(theta);
898 double z = r * sin(theta);
899 return Vector3d(x, y, z);
900 }
901
902
903 /**********************************************************************************
904 * Factory Functions
905 *********************************************************************************/
906
GetRandomIntGenerator(int minval,int maxval,size_t count)907 SeedableIntGeneratorPtr GetRandomIntGenerator(int minval, int maxval, size_t count)
908 {
909 assert (count > 0);
910 Mt19937IntGenerator::ParameterStruct param(minval, maxval);
911 shared_ptr<NumberSequenceFactory<int> > factory = Mt19937IntMetaFactory::GetFactory(param);
912 SeedableIntGeneratorPtr generator(new PrecomputedIntGenerator(factory, count));
913 (void)(*generator)(); // legacy fix
914 return generator;
915 }
916
GetRandomDoubleGenerator(double minval,double maxval,size_t count)917 SeedableDoubleGeneratorPtr GetRandomDoubleGenerator(double minval, double maxval, size_t count)
918 {
919 assert (count > 0);
920 Mt19937DoubleGenerator::ParameterStruct param(minval, maxval);
921 shared_ptr<NumberSequenceFactory<double> > factory(Mt19937DoubleMetaFactory::GetFactory(param));
922 SeedableDoubleGeneratorPtr generator(new PrecomputedDoubleGenerator(factory, count));
923 (void)(*generator)(); // legacy fix
924 return generator;
925 }
926
GetRandomDoubleGenerator(double minval,double maxval)927 SequentialDoubleGeneratorPtr GetRandomDoubleGenerator(double minval, double maxval)
928 {
929 Mt19937DoubleGenerator::ParameterStruct param(minval, maxval);
930 SequentialDoubleGeneratorPtr generator(new Mt19937DoubleGenerator(param));
931 (void)(*generator)(); // legacy fix
932 return generator;
933 }
934
GetIndexedRandomDoubleGenerator(double minval,double maxval,size_t count)935 IndexedDoubleGeneratorPtr GetIndexedRandomDoubleGenerator(double minval, double maxval, size_t count)
936 {
937 assert (count > 0);
938 Mt19937DoubleGenerator::ParameterStruct param(minval, maxval);
939 shared_ptr<NumberSequenceFactory<double> > factory(Mt19937DoubleMetaFactory::GetFactory(param));
940 return IndexedDoubleGeneratorPtr(new PrecomputedDoubleGenerator(factory, count));
941 }
942
GetSubRandomCosWeightedDirectionGenerator(unsigned int id,size_t count)943 SequentialVectorGeneratorPtr GetSubRandomCosWeightedDirectionGenerator(unsigned int id, size_t count)
944 {
945 if ((id == 0) && count && (count < LegacyCosWeightedDirectionGenerator::NumEntries))
946 {
947 LegacyCosWeightedDirectionGenerator::ParameterStruct param;
948 shared_ptr<NumberSequenceFactory<Vector3d> > factory(LegacyCosWeightedDirectionMetaFactory::GetFactory(param));
949 return SequentialVectorGeneratorPtr(new PrecomputedVectorGenerator(factory, count));
950 }
951 else
952 {
953 HaltonCosWeightedDirectionGenerator::ParameterStruct param(primeTable[id % PRIME_TABLE_COUNT], primeTable[(id+1) % PRIME_TABLE_COUNT]);
954 if (count)
955 {
956 shared_ptr<NumberSequenceFactory<Vector3d> > factory(HaltonCosWeightedDirectionMetaFactory::GetFactory(param));
957 return SequentialVectorGeneratorPtr(new PrecomputedVectorGenerator(factory, count));
958 }
959 else
960 return SequentialVectorGeneratorPtr(new HaltonCosWeightedDirectionGenerator(param));
961 }
962 }
963
GetSubRandomDoubleGenerator(unsigned int id,double minval,double maxval,size_t count)964 SequentialDoubleGeneratorPtr GetSubRandomDoubleGenerator(unsigned int id, double minval, double maxval, size_t count)
965 {
966 HaltonDoubleGenerator::ParameterStruct param(primeTable[id % PRIME_TABLE_COUNT], minval, maxval);
967 if (count)
968 {
969 shared_ptr<NumberSequenceFactory<double> > factory(HaltonUniformDoubleMetaFactory::GetFactory(param));
970 return SequentialDoubleGeneratorPtr(new PrecomputedDoubleGenerator(factory, count));
971 }
972 else
973 return SequentialDoubleGeneratorPtr(new HaltonDoubleGenerator(param));
974 }
975
GetSubRandomDirectionGenerator(unsigned int id,size_t count)976 SequentialVectorGeneratorPtr GetSubRandomDirectionGenerator(unsigned int id, size_t count)
977 {
978 HaltonUniformDirectionGenerator::ParameterStruct param(primeTable[id % PRIME_TABLE_COUNT], primeTable[(id+1) % PRIME_TABLE_COUNT]);
979 if (count)
980 {
981 shared_ptr<NumberSequenceFactory<Vector3d> > factory(HaltonUniformDirectionMetaFactory::GetFactory(param));
982 return SequentialVectorGeneratorPtr(new PrecomputedVectorGenerator(factory, count));
983 }
984 else
985 return SequentialVectorGeneratorPtr(new HaltonUniformDirectionGenerator(param));
986 }
987
GetSubRandomOnDiscGenerator(unsigned int id,double radius,size_t count)988 SequentialVector2dGeneratorPtr GetSubRandomOnDiscGenerator(unsigned int id, double radius, size_t count)
989 {
990 HaltonOnDiscGenerator::ParameterStruct param(primeTable[id % PRIME_TABLE_COUNT], primeTable[(id+1) % PRIME_TABLE_COUNT], radius);
991 if (count)
992 {
993 shared_ptr<NumberSequenceFactory<Vector2d> > factory(HaltonOnDiscMetaFactory::GetFactory(param));
994 return SequentialVector2dGeneratorPtr(new PrecomputedVector2dGenerator(factory, count));
995 }
996 else
997 return SequentialVector2dGeneratorPtr(new HaltonOnDiscGenerator(param));
998 }
999
GetSubRandom2dGenerator(unsigned int id,double minX,double maxX,double minY,double maxY,size_t count)1000 SequentialVector2dGeneratorPtr GetSubRandom2dGenerator(unsigned int id, double minX, double maxX, double minY, double maxY, size_t count)
1001 {
1002 Halton2dGenerator::ParameterStruct param(primeTable[id % PRIME_TABLE_COUNT], primeTable[(id+1) % PRIME_TABLE_COUNT], minX, maxX, minY, maxY);
1003 if (count)
1004 {
1005 shared_ptr<NumberSequenceFactory<Vector2d> > factory(Halton2dMetaFactory::GetFactory(param));
1006 return SequentialVector2dGeneratorPtr(new PrecomputedVector2dGenerator(factory, count));
1007 }
1008 else
1009 return SequentialVector2dGeneratorPtr(new Halton2dGenerator(param));
1010 }
1011
1012 } // end of namespace pov
1013