3a9b5dc8 |
1 | // Copyright (c) 2012 Leonhard Gruenschloss (leonhard@gruenschloss.org) |
2 | // |
3 | // Permission is hereby granted, free of charge, to any person obtaining a copy |
4 | // of this software and associated documentation files (the "Software"), to deal |
5 | // in the Software without restriction, including without limitation the rights to |
6 | // use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies |
7 | // of the Software, and to permit persons to whom the Software is furnished to do |
8 | // so, subject to the following conditions: |
9 | // |
10 | // The above copyright notice and this permission notice shall be included in |
11 | // all copies or substantial portions of the Software. |
12 | // |
13 | // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
14 | // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
15 | // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
16 | // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
17 | // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
18 | // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE |
19 | // SOFTWARE. |
20 | |
21 | #ifndef _OpenGl_HaltonSampler_H |
22 | #define _OpenGl_HaltonSampler_H |
23 | |
24 | #include <algorithm> |
25 | #include <vector> |
26 | |
27 | //! Compute points of the Halton sequence with with digit-permutations for different bases. |
28 | class OpenGl_HaltonSampler |
29 | { |
30 | public: |
31 | |
32 | //! Return the number of supported dimensions. |
33 | static unsigned get_num_dimensions() { return 3u; } |
34 | |
35 | public: |
36 | |
37 | //! Init the permutation arrays using Faure-permutations. |
38 | //! Alternatively, initRandom() can be called before the sampling functionality can be used. |
39 | void initFaure(); |
40 | |
41 | //! Init the permutation arrays using randomized permutations. |
42 | //! Alternatively, initFaure() can be called before the sampling functionality can be used. |
43 | //! The client needs to specify a random number generator function object that can be used to generate a random sequence of integers. |
44 | //! That is: if f is a random number generator and N is a positive integer, |
45 | //! then f(N) will return an integer less than N and greater than or equal to 0. |
46 | template <typename Random_number_generator> |
47 | void initRandom (Random_number_generator& theRand); |
48 | |
49 | //! Return the Halton sample for the given dimension (component) and index. |
50 | //! The client must have called initRandom or initFaure() at least once before. |
51 | //! dimension must be smaller than the value returned by get_num_dimensions(). |
52 | float sample (unsigned theDimension, unsigned theIndex) const |
53 | { |
54 | switch (theDimension) |
55 | { |
56 | case 0: return halton2 (theIndex); |
57 | case 1: return halton3 (theIndex); |
58 | case 2: return halton5 (theIndex); |
59 | } |
60 | return 0.0f; |
61 | } |
62 | |
63 | private: |
64 | |
65 | static unsigned short invert (unsigned short theBase, unsigned short theDigits, |
66 | unsigned short theIndex, const std::vector<unsigned short>& thePerm) |
67 | { |
68 | unsigned short aResult = 0; |
69 | for (unsigned short i = 0; i < theDigits; ++i) |
70 | { |
71 | aResult = aResult * theBase + thePerm[theIndex % theBase]; |
72 | theIndex /= theBase; |
73 | } |
74 | return aResult; |
75 | } |
76 | |
77 | void initTables (const std::vector<std::vector<unsigned short> >& thePerm) |
78 | { |
79 | for (unsigned short i = 0; i < 243; ++i) |
80 | { |
81 | myPerm3[i] = invert (3, 5, i, thePerm[3]); |
82 | } |
83 | for (unsigned short i = 0; i < 125; ++i) |
84 | { |
85 | myPerm5[i] = invert (5, 3, i, thePerm[5]); |
86 | } |
87 | } |
88 | |
89 | //! Special case: radical inverse in base 2, with direct bit reversal. |
90 | float halton2 (unsigned theIndex) const |
91 | { |
92 | theIndex = (theIndex << 16) | (theIndex >> 16); |
93 | theIndex = ((theIndex & 0x00ff00ff) << 8) | ((theIndex & 0xff00ff00) >> 8); |
94 | theIndex = ((theIndex & 0x0f0f0f0f) << 4) | ((theIndex & 0xf0f0f0f0) >> 4); |
95 | theIndex = ((theIndex & 0x33333333) << 2) | ((theIndex & 0xcccccccc) >> 2); |
96 | theIndex = ((theIndex & 0x55555555) << 1) | ((theIndex & 0xaaaaaaaa) >> 1); |
97 | union Result { unsigned u; float f; } aResult; // Write reversed bits directly into floating-point mantissa. |
98 | aResult.u = 0x3f800000u | (theIndex >> 9); |
99 | return aResult.f - 1.0f; |
100 | } |
101 | |
102 | float halton3 (unsigned theIndex) const |
103 | { |
104 | return (myPerm3[theIndex % 243u] * 14348907u |
105 | + myPerm3[(theIndex / 243u) % 243u] * 59049u |
106 | + myPerm3[(theIndex / 59049u) % 243u] * 243u |
107 | + myPerm3[(theIndex / 14348907u) % 243u]) * float(0.999999999999999f / 3486784401u); // Results in [0,1). |
108 | } |
109 | |
110 | float halton5 (unsigned theIndex) const |
111 | { |
112 | return (myPerm5[theIndex % 125u] * 1953125u |
113 | + myPerm5[(theIndex / 125u) % 125u] * 15625u |
114 | + myPerm5[(theIndex / 15625u) % 125u] * 125u |
115 | + myPerm5[(theIndex / 1953125u) % 125u]) * float(0.999999999999999f / 244140625u); // Results in [0,1). |
116 | } |
117 | |
118 | private: |
119 | |
120 | unsigned short myPerm3[243]; |
121 | unsigned short myPerm5[125]; |
122 | |
123 | }; |
124 | |
125 | inline void OpenGl_HaltonSampler::initFaure() |
126 | { |
127 | const unsigned THE_MAX_BASE = 5u; |
128 | std::vector<std::vector<unsigned short> > aPerms(THE_MAX_BASE + 1); |
129 | for (unsigned k = 1; k <= 3; ++k) // Keep identity permutations for base 1, 2, 3. |
130 | { |
131 | aPerms[k].resize(k); |
132 | for (unsigned i = 0; i < k; ++i) |
133 | { |
134 | aPerms[k][i] = static_cast<unsigned short> (i); |
135 | } |
136 | } |
137 | |
138 | for (unsigned aBase = 4; aBase <= THE_MAX_BASE; ++aBase) |
139 | { |
140 | aPerms[aBase].resize(aBase); |
141 | const unsigned b = aBase / 2; |
142 | if (aBase & 1) // odd |
143 | { |
144 | for (unsigned i = 0; i < aBase - 1; ++i) |
145 | { |
146 | aPerms[aBase][i + (i >= b)] = aPerms[aBase - 1][i] + (aPerms[aBase - 1][i] >= b); |
147 | } |
148 | aPerms[aBase][b] = static_cast<unsigned short> (b); |
149 | } |
150 | else // even |
151 | { |
152 | for (unsigned i = 0; i < b; ++i) |
153 | { |
154 | aPerms[aBase][i] = 2 * aPerms[b][i]; |
155 | aPerms[aBase][b + i] = 2 * aPerms[b][i] + 1; |
156 | } |
157 | } |
158 | } |
159 | initTables (aPerms); |
160 | } |
161 | |
162 | template <typename Random_number_generator> |
163 | void OpenGl_HaltonSampler::initRandom (Random_number_generator& theRand) |
164 | { |
165 | const unsigned THE_MAX_BASE = 5u; |
166 | std::vector<std::vector<unsigned short> > aPerms(THE_MAX_BASE + 1); |
167 | for (unsigned k = 1; k <= 3; ++k) // Keep identity permutations for base 1, 2, 3. |
168 | { |
169 | aPerms[k].resize (k); |
170 | for (unsigned i = 0; i < k; ++i) |
171 | { |
172 | aPerms[k][i] = i; |
173 | } |
174 | } |
175 | |
176 | for (unsigned aBase = 4; aBase <= THE_MAX_BASE; ++aBase) |
177 | { |
178 | aPerms[aBase].resize (aBase); |
179 | for (unsigned i = 0; i < aBase; ++i) |
180 | { |
181 | aPerms[aBase][i] = i; |
182 | } |
183 | std::random_shuffle (aPerms[aBase].begin(), aPerms[aBase].end(), theRand); |
184 | } |
185 | initTables (aPerms); |
186 | } |
187 | |
188 | #endif // _OpenGl_HaltonSampler_H |