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

RevLine 
[355]1#ifndef __HALTON_H
2#define __HALTON_H
3
[492]4#include <iostream>
5using namespace std;
6
[860]7namespace GtpVisibilityPreprocessor {
8
[1867]9inline float halton(float baseRec, float prev) {
[1876]10  //  float r = 1 - prev - 1e-10f;
11  float r = 1.0f - prev;
[2076]12  //float r = 1.0f - prev;
[1867]13  if (baseRec < r)
14        return prev + baseRec;
[1876]15  float h = baseRec;
16  float hh;
17  do {
18        hh = h;
19        h *= baseRec;
[2076]20  } while (h > r);
[1876]21  return prev + hh + h - 1.0f;
[1867]22}
23
24template<int T>
[2076]25struct Halton {
26  static float _invBases[T];
[1867]27  float _prev[T];
28 
29public:
30 
31  void Reset() {
32        for (int i=0; i < T; i++)
33          _prev[i] = 0;
[492]34  }
[2076]35
36  Halton(const bool initializeBases) {
[1867]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        }
[2076]43  }
44
45  Halton() {
[1867]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
[2076]59struct Halton2 {
[1867]60  static float _invBases[2];
61  float _prev[2];
62 
[492]63public:
64 
65  void Reset() {
66        _prev[0] =_prev[1] = 0;
67  }
68 
69  Halton2() {
[497]70        _invBases[0] = 1.0f/2;
71        _invBases[1] = 1.0f/3;
[492]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};
[355]83
[492]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;
[355]95        }
[492]96  }
97        if(number == 2) {
98          isIt = false;
99        }
100        return isIt;
101}
102
[2076]103 
104  /**
[492]105 * Find the nth prime number.
106 * @param index the ordinal position in the sequence
107 */
108inline int FindPrime(const int index) {
[1867]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
[492]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 {
[355]130public:
[492]131  int index;
[355]132
[2105]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
[492]141  HaltonSequence():index(1) {}
[355]142
[492]143  void Reset() {
144        index = 1;
145  }
146
[1876]147  void
[2076]148  GetNext(const int dimensions, float *p);
[1876]149 
[492]150  void GenerateNext() {
151        index++;
152  }
153 
[2076]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 
[492]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
[2076]179  double GetNumberOld(const int dimension)  {
[492]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;
[497]189        if((base >= 2) && (index >= 1)) {
[492]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;
[355]198        }
[492]199        else {
200          cerr<<"Error generating Halton sequence."<<endl;
201          exit(1);
[355]202        }
[492]203  }
[355]204};
205
206extern Halton2 halton2;
[860]207}
[355]208
209#endif
Note: See TracBrowser for help on using the repository browser.