1 // Copyright (c) 2012 Leonhard Gruenschloss (leonhard@gruenschloss.org)
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:
10 // The above copyright notice and this permission notice shall be included in
11 // all copies or substantial portions of the Software.
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
21 #ifndef _OpenGl_HaltonSampler_H
22 #define _OpenGl_HaltonSampler_H
27 //! Compute points of the Halton sequence with digit-permutations for different bases.
28 class OpenGl_HaltonSampler
32 //! Return the number of supported dimensions.
33 static unsigned get_num_dimensions() { return 3u; }
37 //! Init the permutation arrays using Faure-permutations.
38 OpenGl_HaltonSampler()
43 //! Return the Halton sample for the given dimension (component) and index.
44 //! The client must have called initFaure() at least once before.
45 //! dimension must be smaller than the value returned by get_num_dimensions().
46 float sample (unsigned theDimension, unsigned theIndex) const
50 case 0: return halton2 (theIndex);
51 case 1: return halton3 (theIndex);
52 case 2: return halton5 (theIndex);
59 //! Init the permutation arrays using Faure-permutations.
62 static unsigned short invert (unsigned short theBase, unsigned short theDigits,
63 unsigned short theIndex, const std::vector<unsigned short>& thePerm)
65 unsigned short aResult = 0;
66 for (unsigned short i = 0; i < theDigits; ++i)
68 aResult = aResult * theBase + thePerm[theIndex % theBase];
74 void initTables (const std::vector<std::vector<unsigned short> >& thePerm)
76 for (unsigned short i = 0; i < 243; ++i)
78 myPerm3[i] = invert (3, 5, i, thePerm[3]);
80 for (unsigned short i = 0; i < 125; ++i)
82 myPerm5[i] = invert (5, 3, i, thePerm[5]);
86 //! Special case: radical inverse in base 2, with direct bit reversal.
87 float halton2 (unsigned theIndex) const
89 theIndex = (theIndex << 16) | (theIndex >> 16);
90 theIndex = ((theIndex & 0x00ff00ff) << 8) | ((theIndex & 0xff00ff00) >> 8);
91 theIndex = ((theIndex & 0x0f0f0f0f) << 4) | ((theIndex & 0xf0f0f0f0) >> 4);
92 theIndex = ((theIndex & 0x33333333) << 2) | ((theIndex & 0xcccccccc) >> 2);
93 theIndex = ((theIndex & 0x55555555) << 1) | ((theIndex & 0xaaaaaaaa) >> 1);
94 union Result { unsigned u; float f; } aResult; // Write reversed bits directly into floating-point mantissa.
95 aResult.u = 0x3f800000u | (theIndex >> 9);
96 return aResult.f - 1.0f;
99 float halton3 (unsigned theIndex) const
101 return (myPerm3[theIndex % 243u] * 14348907u
102 + myPerm3[(theIndex / 243u) % 243u] * 59049u
103 + myPerm3[(theIndex / 59049u) % 243u] * 243u
104 + myPerm3[(theIndex / 14348907u) % 243u]) * float(0.999999999999999 / 3486784401u); // Results in [0,1).
107 float halton5 (unsigned theIndex) const
109 return (myPerm5[theIndex % 125u] * 1953125u
110 + myPerm5[(theIndex / 125u) % 125u] * 15625u
111 + myPerm5[(theIndex / 15625u) % 125u] * 125u
112 + myPerm5[(theIndex / 1953125u) % 125u]) * float(0.999999999999999 / 244140625u); // Results in [0,1).
117 unsigned short myPerm3[243];
118 unsigned short myPerm5[125];
122 inline void OpenGl_HaltonSampler::initFaure()
124 const unsigned THE_MAX_BASE = 5u;
125 std::vector<std::vector<unsigned short> > aPerms(THE_MAX_BASE + 1);
126 for (unsigned k = 1; k <= 3; ++k) // Keep identity permutations for base 1, 2, 3.
129 for (unsigned i = 0; i < k; ++i)
131 aPerms[k][i] = static_cast<unsigned short> (i);
135 for (unsigned aBase = 4; aBase <= THE_MAX_BASE; ++aBase)
137 aPerms[aBase].resize(aBase);
138 const unsigned b = aBase / 2;
139 if (aBase & 1) // odd
141 for (unsigned i = 0; i < aBase - 1; ++i)
143 aPerms[aBase][i + (i >= b)] = aPerms[aBase - 1][i] + (aPerms[aBase - 1][i] >= b);
145 aPerms[aBase][b] = static_cast<unsigned short> (b);
149 for (unsigned i = 0; i < b; ++i)
151 aPerms[aBase][i] = 2 * aPerms[b][i];
152 aPerms[aBase][b + i] = 2 * aPerms[b][i] + 1;
159 #endif // _OpenGl_HaltonSampler_H