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

Revision 2176, 3.8 KB checked in by mattausch, 18 years ago (diff)

removed using namespace std from .h

Line 
1#ifndef __HALTON_H
2#define __HALTON_H
3
4#include <iostream>
5
6namespace GtpVisibilityPreprocessor {
7
8inline float halton(float baseRec, float prev) {
9  //  float r = 1 - prev - 1e-10f;
10  float r = 1.0f - prev;
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>
24struct Halton {
25  static 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(const bool initializeBases) {
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  }
43
44  Halton() {
45        Reset();
46  }
47 
48  void
49  GetNext(float *a) {
50        for (int i=0; i < T; i++) {
51          a[i] = halton(_invBases[i], _prev[i]);
52          _prev[i] = a[i];
53        }
54  }
55 
56};
57
58struct Halton2 {
59  static float _invBases[2];
60  float _prev[2];
61 
62public:
63 
64  void Reset() {
65        _prev[0] =_prev[1] = 0;
66  }
67 
68  Halton2() {
69        _invBases[0] = 1.0f/2;
70        _invBases[1] = 1.0f/3;
71        Reset();
72  }
73 
74  void
75  GetNext(float &a, float &b) {
76        a = halton(_invBases[0], _prev[0]);
77        b = halton(_invBases[1], _prev[1]);
78        _prev[0] = a;
79        _prev[1] = b;
80  }
81};
82
83
84/**
85 * Assert whether the argument is a prime number.
86 * @param number the number to be checked
87 */
88inline bool IsPrime(const int number) {
89  bool isIt = true;
90  for(int i = 2; i < number; i++) {
91        if(number % i == 0) {
92          isIt = false;
93          break;
94        }
95  }
96        if(number == 2) {
97          isIt = false;
98        }
99        return isIt;
100}
101
102 
103  /**
104 * Find the nth prime number.
105 * @param index the ordinal position in the sequence
106 */
107inline int FindPrime(const int index) {
108  //  if (index < 1) {
109  //    cerr<<"FindPrime: The argument must be non-negative."<<endl;
110  //    return -1;
111  //  }
112
113  const int primes[] = {-1, 1, 3, 5, 7, 11, 13};
114  if (index <= 6)
115        return primes[index];
116
117  int prime = 1;
118  int found = 1;
119  while(found != index) {
120        prime += 2;
121          if(IsPrime(prime) == true) {
122                found++;
123          }
124  }
125  return prime;
126}
127
128struct HaltonSequence {
129public:
130  int index;
131
132  static int sPregeneratedDim;
133  static int sPregeneratedNumber;
134  static float *sPregeneratedValues;
135 
136  // special construtor for pregenerating static halton sequences
137  HaltonSequence(const int dim,
138                                 const int number);
139
140  HaltonSequence():index(1) {}
141
142  void Reset() {
143        index = 1;
144  }
145
146  void
147  GetNext(const int dimensions, float *p);
148 
149  void GenerateNext() {
150        index++;
151  }
152 
153  double GetNumber(const int dimension)  {
154        int base = FindPrime(dimension);
155        if(base == 1) {
156          base++;  //The first dimension uses base 2.
157        }
158       
159        int _p1 = base;
160        float _ip1 = 1.0f/base;
161        float p, u=0.0f;
162    int kk, a;
163       
164    // the first coordinate
165    for (p = _ip1, kk = index ; kk ;  p *= _ip1, kk /= _p1)   
166      if ((a = kk % _p1))
167                u += a * p;
168       
169        return u;
170  }
171 
172  /**
173   * Returns the nth number in the sequence, taken from a specified dimension.
174   * @param index the ordinal position in the sequence
175   * @param dimension the dimension
176   */
177
178  double GetNumberOld(const int dimension)  {
179        int base = FindPrime(dimension);
180        if(base == 1) {
181          base++;  //The first dimension uses base 2.
182        }
183        double remainder;
184        double output = 0.0;
185        double fraction = 1.0 / (double)base;
186        int N1 = 0;
187        int copyOfIndex = index;
188        if((base >= 2) && (index >= 1)) {
189          while(copyOfIndex > 0) {
190                N1 = (copyOfIndex / base);
191                remainder = copyOfIndex % base;
192                output += fraction * remainder;
193                copyOfIndex = (int)(copyOfIndex / base);
194                fraction /= (double)base;
195          }
196          return output;
197        }
198        else {
199          std::cerr<<"Error generating Halton sequence."<<std::endl;
200          exit(1);
201        }
202  }
203};
204
205extern Halton2 halton2;
206}
207
208#endif
Note: See TracBrowser for help on using the repository browser.