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

Revision 2105, 3.8 KB checked in by bittner, 18 years ago (diff)

simple ray separated

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