Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-14 09:15:37

0001 // Copyright (c) 2012 Leonhard Gruenschloss (leonhard@gruenschloss.org)
0002 //
0003 // Permission is hereby granted, free of charge, to any person obtaining a copy
0004 // of this software and associated documentation files (the "Software"), to deal
0005 // in the Software without restriction, including without limitation the rights to
0006 // use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
0007 // of the Software, and to permit persons to whom the Software is furnished to do
0008 // so, subject to the following conditions:
0009 //
0010 // The above copyright notice and this permission notice shall be included in
0011 // all copies or substantial portions of the Software.
0012 //
0013 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
0014 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0015 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
0016 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0017 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
0018 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
0019 // SOFTWARE.
0020 
0021 #ifndef _OpenGl_HaltonSampler_H
0022 #define _OpenGl_HaltonSampler_H
0023 
0024 #include <NCollection_LinearVector.hxx>
0025 
0026 //! Compute points of the Halton sequence with digit-permutations for different bases.
0027 class OpenGl_HaltonSampler
0028 {
0029 public:
0030   //! Return the number of supported dimensions.
0031   static unsigned get_num_dimensions() { return 3u; }
0032 
0033 public:
0034   //! Init the permutation arrays using Faure-permutations.
0035   OpenGl_HaltonSampler() { initFaure(); }
0036 
0037   //! Return the Halton sample for the given dimension (component) and index.
0038   //! The client must have called initFaure() at least once before.
0039   //! dimension must be smaller than the value returned by get_num_dimensions().
0040   float sample(unsigned theDimension, unsigned theIndex) const
0041   {
0042     switch (theDimension)
0043     {
0044       case 0:
0045         return halton2(theIndex);
0046       case 1:
0047         return halton3(theIndex);
0048       case 2:
0049         return halton5(theIndex);
0050     }
0051     return 0.0f;
0052   }
0053 
0054 private:
0055   //! Init the permutation arrays using Faure-permutations.
0056   void initFaure();
0057 
0058   static unsigned short invert(unsigned short                                  theBase,
0059                                unsigned short                                  theDigits,
0060                                unsigned short                                  theIndex,
0061                                const NCollection_LinearVector<unsigned short>& thePerm)
0062   {
0063     unsigned short aResult = 0;
0064     for (unsigned short i = 0; i < theDigits; ++i)
0065     {
0066       aResult = aResult * theBase + thePerm[theIndex % theBase];
0067       theIndex /= theBase;
0068     }
0069     return aResult;
0070   }
0071 
0072   void initTables(const NCollection_LinearVector<NCollection_LinearVector<unsigned short>>& thePerm)
0073   {
0074     for (unsigned short i = 0; i < 243; ++i)
0075     {
0076       myPerm3[i] = invert(3, 5, i, thePerm[3]);
0077     }
0078     for (unsigned short i = 0; i < 125; ++i)
0079     {
0080       myPerm5[i] = invert(5, 3, i, thePerm[5]);
0081     }
0082   }
0083 
0084   //! Special case: radical inverse in base 2, with direct bit reversal.
0085   float halton2(unsigned theIndex) const
0086   {
0087     theIndex = (theIndex << 16) | (theIndex >> 16);
0088     theIndex = ((theIndex & 0x00ff00ff) << 8) | ((theIndex & 0xff00ff00) >> 8);
0089     theIndex = ((theIndex & 0x0f0f0f0f) << 4) | ((theIndex & 0xf0f0f0f0) >> 4);
0090     theIndex = ((theIndex & 0x33333333) << 2) | ((theIndex & 0xcccccccc) >> 2);
0091     theIndex = ((theIndex & 0x55555555) << 1) | ((theIndex & 0xaaaaaaaa) >> 1);
0092 
0093     // clang-format off
0094     union Result { unsigned u; float f; } aResult; // Write reversed bits directly into floating-point mantissa.
0095 
0096     // clang-format on
0097     aResult.u = 0x3f800000u | (theIndex >> 9);
0098     return aResult.f - 1.0f;
0099   }
0100 
0101   float halton3(unsigned theIndex) const
0102   {
0103     return (myPerm3[theIndex % 243u] * 14348907u + myPerm3[(theIndex / 243u) % 243u] * 59049u
0104             + myPerm3[(theIndex / 59049u) % 243u] * 243u + myPerm3[(theIndex / 14348907u) % 243u])
0105            * float(0.999999999999999 / 3486784401u); // Results in [0,1).
0106   }
0107 
0108   float halton5(unsigned theIndex) const
0109   {
0110     return (myPerm5[theIndex % 125u] * 1953125u + myPerm5[(theIndex / 125u) % 125u] * 15625u
0111             + myPerm5[(theIndex / 15625u) % 125u] * 125u + myPerm5[(theIndex / 1953125u) % 125u])
0112            * float(0.999999999999999 / 244140625u); // Results in [0,1).
0113   }
0114 
0115 private:
0116   unsigned short myPerm3[243];
0117   unsigned short myPerm5[125];
0118 };
0119 
0120 inline void OpenGl_HaltonSampler::initFaure()
0121 {
0122   const unsigned                                                     THE_MAX_BASE = 5u;
0123   NCollection_LinearVector<NCollection_LinearVector<unsigned short>> aPerms;
0124   aPerms.Resize(THE_MAX_BASE + 1);
0125   for (unsigned k = 1; k <= 3; ++k) // Keep identity permutations for base 1, 2, 3.
0126   {
0127     aPerms[k].Resize(k);
0128     for (unsigned i = 0; i < k; ++i)
0129     {
0130       aPerms[k][i] = static_cast<unsigned short>(i);
0131     }
0132   }
0133 
0134   for (unsigned aBase = 4; aBase <= THE_MAX_BASE; ++aBase)
0135   {
0136     aPerms[aBase].Resize(aBase);
0137     const unsigned b = aBase / 2;
0138     if (aBase & 1) // odd
0139     {
0140       for (unsigned i = 0; i < aBase - 1; ++i)
0141       {
0142         aPerms[aBase][i + (i >= b)] = aPerms[aBase - 1][i] + (aPerms[aBase - 1][i] >= b);
0143       }
0144       aPerms[aBase][b] = static_cast<unsigned short>(b);
0145     }
0146     else // even
0147     {
0148       for (unsigned i = 0; i < b; ++i)
0149       {
0150         aPerms[aBase][i]     = 2 * aPerms[b][i];
0151         aPerms[aBase][b + i] = 2 * aPerms[b][i] + 1;
0152       }
0153     }
0154   }
0155   initTables(aPerms);
0156 }
0157 
0158 #endif // _OpenGl_HaltonSampler_H