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

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