source: GTP/trunk/App/Demos/Vis/FriendlyCulling/src/Halton.h @ 3129

#ifndef __HALTON_H
#define __HALTON_H

#include <iostream>


/** Assert whether the argument is a prime number.
        @param number the number to be checked
*/
inline bool IsPrime(const int number) 
{
        bool isIt = true;

        for(int i = 2; i < number; ++ i) 
        {
                if(number % i == 0) 
                {
                        isIt = false;
                        break;
                }
        }

        if(number == 2) 
        {
                isIt = false;
        }

        return isIt;
}


/**
        Find the nth prime number.
        @param index the ordinal position in the sequence
*/
inline int FindPrime(const int index) 
{

        const int primes[] = {-1, 1, 3, 5, 7, 11, 13};
        if (index <= 6)
                return primes[index];

        int prime = 1;
        int found = 1;
        while(found != index) {
                prime += 2;
                if(IsPrime(prime) == true) {
                        found++;
                }
        }
        return prime;
}

  
inline float halton(float baseRec, float prev) 
{
        float r = 1.0f - prev;
        
        if (baseRec < r)
                return prev + baseRec;
        
        float h = baseRec;
        
        float hh;
        
        do 
        {
                hh = h;
                h *= baseRec;
        } while (h > r);
        
        return prev + hh + h - 1.0f;
}


template<int T> struct Halton 
{
        static float _invBases[T];
        float _prev[T];

public:

        void Reset() 
        {
                for (int i=0; i < T; i++) 
                        _prev[i] = 0;
        }

        Halton(const bool initializeBases) 
        {
                for (int i=0; i < T; i++) 
                {
                        int base = FindPrime(i+1);
                
                        if (base == 1)
                                base++;
                        _invBases[i] = 1.0f/base;
                }
        }

        Halton() 
        {
                Reset();
        }
  
        void GetNext(float *a) 
        {
                for (int i=0; i < T; i++) 
                {
                        a[i] = halton(_invBases[i], _prev[i]);
                        _prev[i] = a[i];
                }
        }

};

struct Halton2 
{
        static float _invBases[2];
        float _prev[2];

public:

        void Reset() 
        {
                _prev[0] =_prev[1] = 0;
        }

        Halton2() 
        {
                _invBases[0] = 1.0f / 2;
                _invBases[1] = 1.0f / 3;
                Reset();
        }

        void GetNext(float &a, float &b) 
        {
                a = halton(_invBases[0], _prev[0]);
                b = halton(_invBases[1], _prev[1]);

                _prev[0] = a;
                _prev[1] = b;
        }
};


  
struct HaltonSequence 
{
public:
        int index;

        static int sPregeneratedDim;
        static int sPregeneratedNumber;
        static float *sPregeneratedValues;

        // special construtor for pregenerating static halton sequences
        HaltonSequence(int dim, int number);

        HaltonSequence(): index(1) 
        { 
                float dummy[2]; 
                for (int i = 0; i < 20; ++ i) GetNext(2, dummy);
        }

        void Reset() { index = 1; }

        void GetNext(int dimensions, float *p);

        void GenerateNext() { ++ index; }

        double GetNumber(const int dimension)  
        {
                int base = FindPrime(dimension);
                if (base == 1) 
                {
                        base++;  // The first dimension uses base 2.
                } 

                int _p1 = base;
                float _ip1 = 1.0f / base;
                float p, u=0.0f;
                int kk, a;

                // the first coordinate
                for (p = _ip1, kk = index ; kk ;  p *= _ip1, kk /= _p1)   
                        if ((a = kk % _p1))
                                u += a * p;

                return u;
        }

        /**
                Returns the nth number in the sequence, taken from a specified dimension.
                @param index the ordinal position in the sequence
                @param dimension the dimension
        */
        double GetNumberOld(int dimension)  
        {
                int base = FindPrime(dimension);
                if(base == 1) 
                {
                        ++ base;  //The first dimension uses base 2.
                } 

                double remainder;
                double output = 0.0;
                double fraction = 1.0 / (double)base;
                int N1 = 0;
                int copyOfIndex = index;
                
                if ((base >= 2) && (index >= 1)) 
                {
                        while (copyOfIndex > 0) 
                        {
                                N1 = (copyOfIndex / base);
                                remainder = copyOfIndex % base;
                                output += fraction * remainder;
                                copyOfIndex = (int)(copyOfIndex / base);
                                fraction /= (double)base;
                        } 
                        return output;
                }
                else 
                {
                        std::cerr << "Error generating Halton sequence." << std::endl;
                        exit(1);
                }
        }
};


#endif