source: GTP/trunk/Lib/Vis/Preprocessing/src/Halton.h @ 2575

#ifndef __HALTON_H
#define __HALTON_H

#include <iostream>

namespace GtpVisibilityPreprocessor {


/**
 * 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) {
  //  if (index < 1) {
  //    cerr<<"FindPrime: The argument must be non-negative."<<endl;
  //    return -1;
  //  }

  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 - prev - 1e-10f;
  float r = 1.0f - 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(const int dim,
                                 const int number);

  HaltonSequence():index(1) {}

  void Reset() {
        index = 1;
  }

  void
  GetNext(const 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(const 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);
        }
  }
};

extern Halton2 halton2;
}

#endif