Changeset 492 for trunk/VUT/GtpVisibilityPreprocessor/src/Halton.h

Timestamp:

01/03/06 23:33:45 (18 years ago)

Author:

bittner

Message:

Large merge - viewcells seem not functional now

File:

• trunk/VUT/GtpVisibilityPreprocessor/src/Halton.h (modified) (1 diff)

trunk/VUT/GtpVisibilityPreprocessor/src/Halton.h

@@ -2 +2 @@
 #define __HALTON_H
 
+#include <iostream>
+using namespace std;
+
 class Halton2 {
-        
-        float _invBases[2];
-        float _prev[2];
+  static float _invBases[2];
+  float _prev[2];
+  
+  float halton(float baseRec, float prev) const {
+        float r = 1 - prev - 1e-10;
+        if (baseRec < r)
+          return prev + baseRec;
+        float h = baseRec;
+        float hh;
+        do {
+          hh = h;
+          h *= baseRec;
+        } while (h >= r);
+        return prev + hh + h - 1;
+  }
+  
+public:
+  
+  void Reset() {
+        _prev[0] =_prev[1] = 0;
+  }
+  
+  Halton2() {
+        _invBases[0] = 1./2;
+        _invBases[1] = 1./3;
+        Reset();
+  }
+  
+  void
+  GetNext(float &a, float &b) {
+        a = halton(_invBases[0], _prev[0]);
+        b = halton(_invBases[1], _prev[1]);
+        _prev[0] = a;
+        _prev[1] = b;
+  }
+};
 
-        float halton(float baseRec, float prev) const {
-                float r = 1 - prev - 1e-10;
-                if (baseRec < r) return prev + baseRec;
-                float h = baseRec;
-                float hh;
-                do {
-                        hh = h;
-                        h *= baseRec;
-                } while (h >= r);
-                return prev + hh + h - 1;
+
+/**
+ * Assert whether the argument is a prime number.
+ * @param number the number to be checked
+ */
+inline bool IsPrime(const int number) {
+  bool isIt = true;
+  for(int i = 2; i < number; i++) {
+        if(number % i == 0) {
+          isIt = false;
+          break;
         }
-        
+  }
+        if(number == 2) {
+          isIt = false;
+        }
+        return isIt;
+}
+
+/**
+ * Find the nth prime number.
+ * @param index the ordinal position in the sequence
+ */
+inline int FindPrime(const int index) {
+  if(index < 1) {
+        cerr<<"FindPrime: The argument must be non-negative."<<endl;
+        return -1;
+  }
+  int prime = 1;
+  int found = 1;
+  while(found != index) {
+        prime += 2;
+          if(IsPrime(prime) == true) {
+                found++;
+          }
+  }
+  return prime;
+}
+
+struct HaltonSequence {
 public:
+  int index;
 
-        void Reset() {
-                _prev[0] =_prev[1] = 0;
+  HaltonSequence():index(1) {}
+
+  void Reset() {
+        index = 1;
+  }
+
+  void GenerateNext() {
+        index++;
+  }
+  
+  /**
+   * Returns the nth number in the sequence, taken from a specified dimension.
+   * @param index the ordinal position in the sequence
+   * @param dimension the dimension
+   */
+
+  double GetNumber(const int dimension)  {
+        int base = FindPrime(dimension);
+        if(base == 1) {
+          base++;  //The first dimension uses base 2.
+        } 
+        double remainder;
+        double output = 0.0;
+        double fraction = 1.0 / (double)base;
+        int N1 = 0;
+        int copyOfIndex = index;
+        if(base >= 2 & index >= 1) {
+          while(copyOfIndex > 0) {
+                N1 = (copyOfIndex / base);
+                remainder = copyOfIndex % base;
+                output += fraction * remainder;
+                copyOfIndex = (int)(copyOfIndex / base);
+                fraction /= (double)base;
+          } 
+          return output;
         }
-
-        Halton2() {
-                _invBases[0] = 1./2;
-                _invBases[1] = 1./3;
-                Reset();
+        else {
+          cerr<<"Error generating Halton sequence."<<endl;
+          exit(1);
         }
-        
-        void
-        GetNext(float &a, float &b) {
-                a = halton(_invBases[0], _prev[0]);
-                b = halton(_invBases[1], _prev[1]);
-                _prev[0] = a;
-                _prev[1] = b;
-        }
+  }
 };