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

Revision 497, 2.5 KB checked in by mattausch, 19 years ago (diff)

beware: bug in view cells merging

Line 
1#ifndef __HALTON_H
2#define __HALTON_H
3
4#include <iostream>
5using namespace std;
6
7class Halton2 {
8  static float _invBases[2];
9  float _prev[2];
10 
11  float halton(float baseRec, float prev) const {
12        float r = 1 - prev - 1e-10f;
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;
22  }
23 
24public:
25 
26  void Reset() {
27        _prev[0] =_prev[1] = 0;
28  }
29 
30  Halton2() {
31        _invBases[0] = 1.0f/2;
32        _invBases[1] = 1.0f/3;
33        Reset();
34  }
35 
36  void
37  GetNext(float &a, float &b) {
38        a = halton(_invBases[0], _prev[0]);
39        b = halton(_invBases[1], _prev[1]);
40        _prev[0] = a;
41        _prev[1] = b;
42  }
43};
44
45
46/**
47 * Assert whether the argument is a prime number.
48 * @param number the number to be checked
49 */
50inline bool IsPrime(const int number) {
51  bool isIt = true;
52  for(int i = 2; i < number; i++) {
53        if(number % i == 0) {
54          isIt = false;
55          break;
56        }
57  }
58        if(number == 2) {
59          isIt = false;
60        }
61        return isIt;
62}
63
64/**
65 * Find the nth prime number.
66 * @param index the ordinal position in the sequence
67 */
68inline int FindPrime(const int index) {
69  if(index < 1) {
70        cerr<<"FindPrime: The argument must be non-negative."<<endl;
71        return -1;
72  }
73  int prime = 1;
74  int found = 1;
75  while(found != index) {
76        prime += 2;
77          if(IsPrime(prime) == true) {
78                found++;
79          }
80  }
81  return prime;
82}
83
84struct HaltonSequence {
85public:
86  int index;
87
88  HaltonSequence():index(1) {}
89
90  void Reset() {
91        index = 1;
92  }
93
94  void GenerateNext() {
95        index++;
96  }
97 
98  /**
99   * Returns the nth number in the sequence, taken from a specified dimension.
100   * @param index the ordinal position in the sequence
101   * @param dimension the dimension
102   */
103
104  double GetNumber(const int dimension)  {
105        int base = FindPrime(dimension);
106        if(base == 1) {
107          base++;  //The first dimension uses base 2.
108        }
109        double remainder;
110        double output = 0.0;
111        double fraction = 1.0 / (double)base;
112        int N1 = 0;
113        int copyOfIndex = index;
114        if((base >= 2) && (index >= 1)) {
115        // changed by matt
116        //if(base >= 2 & index >= 1) {
117          while(copyOfIndex > 0) {
118                N1 = (copyOfIndex / base);
119                remainder = copyOfIndex % base;
120                output += fraction * remainder;
121                copyOfIndex = (int)(copyOfIndex / base);
122                fraction /= (double)base;
123          }
124          return output;
125        }
126        else {
127          cerr<<"Error generating Halton sequence."<<endl;
128          exit(1);
129        }
130  }
131};
132
133extern Halton2 halton2;
134
135#endif
Note: See TracBrowser for help on using the repository browser.