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

Revision 860, 2.5 KB checked in by mattausch, 19 years ago (diff)
Line 
1#ifndef __HALTON_H
2#define __HALTON_H
3
4#include <iostream>
5using namespace std;
6
7namespace GtpVisibilityPreprocessor {
8
9class Halton2 {
10  static float _invBases[2];
11  float _prev[2];
12 
13  float halton(float baseRec, float prev) const {
14        float r = 1 - prev - 1e-10f;
15        if (baseRec < r)
16          return prev + baseRec;
17        float h = baseRec;
18        float hh;
19        do {
20          hh = h;
21          h *= baseRec;
22        } while (h >= r);
23        return prev + hh + h - 1;
24  }
25 
26public:
27 
28  void Reset() {
29        _prev[0] =_prev[1] = 0;
30  }
31 
32  Halton2() {
33        _invBases[0] = 1.0f/2;
34        _invBases[1] = 1.0f/3;
35        Reset();
36  }
37 
38  void
39  GetNext(float &a, float &b) {
40        a = halton(_invBases[0], _prev[0]);
41        b = halton(_invBases[1], _prev[1]);
42        _prev[0] = a;
43        _prev[1] = b;
44  }
45};
46
47
48/**
49 * Assert whether the argument is a prime number.
50 * @param number the number to be checked
51 */
52inline bool IsPrime(const int number) {
53  bool isIt = true;
54  for(int i = 2; i < number; i++) {
55        if(number % i == 0) {
56          isIt = false;
57          break;
58        }
59  }
60        if(number == 2) {
61          isIt = false;
62        }
63        return isIt;
64}
65
66/**
67 * Find the nth prime number.
68 * @param index the ordinal position in the sequence
69 */
70inline int FindPrime(const int index) {
71  if(index < 1) {
72        cerr<<"FindPrime: The argument must be non-negative."<<endl;
73        return -1;
74  }
75  int prime = 1;
76  int found = 1;
77  while(found != index) {
78        prime += 2;
79          if(IsPrime(prime) == true) {
80                found++;
81          }
82  }
83  return prime;
84}
85
86struct HaltonSequence {
87public:
88  int index;
89
90  HaltonSequence():index(1) {}
91
92  void Reset() {
93        index = 1;
94  }
95
96  void GenerateNext() {
97        index++;
98  }
99 
100  /**
101   * Returns the nth number in the sequence, taken from a specified dimension.
102   * @param index the ordinal position in the sequence
103   * @param dimension the dimension
104   */
105
106  double GetNumber(const int dimension)  {
107        int base = FindPrime(dimension);
108        if(base == 1) {
109          base++;  //The first dimension uses base 2.
110        }
111        double remainder;
112        double output = 0.0;
113        double fraction = 1.0 / (double)base;
114        int N1 = 0;
115        int copyOfIndex = index;
116        if((base >= 2) && (index >= 1)) {
117        // changed by matt
118        //if(base >= 2 & index >= 1) {
119          while(copyOfIndex > 0) {
120                N1 = (copyOfIndex / base);
121                remainder = copyOfIndex % base;
122                output += fraction * remainder;
123                copyOfIndex = (int)(copyOfIndex / base);
124                fraction /= (double)base;
125          }
126          return output;
127        }
128        else {
129          cerr<<"Error generating Halton sequence."<<endl;
130          exit(1);
131        }
132  }
133};
134
135extern Halton2 halton2;
136}
137
138#endif
Note: See TracBrowser for help on using the repository browser.