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