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

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