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

Revision 2575, 3.8 KB checked in by bittner, 17 years ago (diff)

big merge: preparation for havran ray caster, check if everything works

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