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

#ifndef __HALTON_H
#define __HALTON_H

#include <iostream>
using namespace std;

namespace GtpVisibilityPreprocessor {

class Halton2 {
  static float _invBases[2];
  float _prev[2];
  
  float halton(float baseRec, float prev) const {
        float r = 1 - prev - 1e-10f;
        if (baseRec < r)
          return prev + baseRec;
        float h = baseRec;
        float hh;
        do {
          hh = h;
          h *= baseRec;
        } while (h >= r);
        return prev + hh + h - 1;
  }
  
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;
  }
};


/**
 * 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;
  }
  int prime = 1;
  int found = 1;
  while(found != index) {
        prime += 2;
          if(IsPrime(prime) == true) {
                found++;
          }
  }
  return prime;
}

struct HaltonSequence {
public:
  int index;

  HaltonSequence():index(1) {}

  void Reset() {
        index = 1;
  }

  void GenerateNext() {
        index++;
  }
  
  /**
   * 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 GetNumber(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)) {
        // changed by matt
        //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 {
          cerr<<"Error generating Halton sequence."<<endl;
          exit(1);
        }
  }
};

extern Halton2 halton2;
}

#endif