0027607: Visualization - Implement adaptive screen space sampling in path tracing
[occt.git] / src / OpenGl / OpenGl_HaltonSampler.hxx
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