improved normalization

gravitino · Mar 13, 2015 · 35fce57 · 35fce57
1 parent 320ad7d
commit 35fce57
Show file tree

Hide file tree

Showing 3 changed files with 19 additions and 16 deletions.
diff --git a/blending.cpp b/blending.cpp
@@ -3,22 +3,22 @@
 
 int main () {
 
-	using namespace average;
+    using namespace average;
 
     std::mt19937 engine;
     engine.seed(0);
     std::uniform_real_distribution<double> dist(-1, 1);
 
     size_t N = 10000;
 
-	auto q0 = quat<double>(0, 1, 1, 1).exp(); q0.print();
+    auto q0 = quat<double>(0, 1, 1, 1).exp(); q0.print();
 
     std::vector<double> weights(N, 1.0/N);
     std::vector<quat<double>> quats;
 
     for (size_t i = 0; i < N; i++)
-    	quats.push_back( q0*quat<double>(1+dist(engine), dist(engine), 
-    		                               dist(engine), dist(engine)).N());
+        quats.push_back( q0*quat<double>(1+dist(engine), dist(engine), 
+                                           dist(engine), dist(engine)).N());
 
     QLA(quats).print();
     std::cout << QLA(quats).logdist(q0) << std::endl;
@@ -36,10 +36,10 @@ int main () {
     std::vector<dualquat<double>> Quats;
 
     for (size_t i = 0; i < N; i++)
-    	Quats.push_back(Q0*dualquat<double>(1+dist(engine), dist(engine), 
-    		                                  dist(engine), dist(engine),
-    		                                  dist(engine), dist(engine),
-    		                                  dist(engine), dist(engine)).N());
+        Quats.push_back(Q0*dualquat<double>(1+dist(engine), dist(engine), 
+                                              dist(engine), dist(engine),
+                                              dist(engine), dist(engine),
+                                              dist(engine), dist(engine)).N());
 
     DLA(Quats).print();
     std::cout << DLA(Quats).logdist(Q0) << std::endl;

diff --git a/dualquat.hpp b/dualquat.hpp
@@ -89,7 +89,7 @@ struct quat {
 
         assert(w*w < EPS(value_t));
 
-        value_t dst = std::sqrt(x*x+y*y+z*z+EPS(value_t));
+        value_t dst = std::sqrt(x*x+y*y+z*z);
         value_t fac = std::sin(dst)/dst;        
 
         return quat<value_t>(std::cos(dst), x*fac, y*fac, z*fac);
@@ -113,8 +113,11 @@ struct quat {
 
         assert(isunit());
 
-        const value_t inv = 1.0/std::sqrt(x*x+y*y+z*z+EPS(value_t));
-        const value_t fac = std::acos(w > 1 ? 1 : w < -1 ? -1 : w)*inv;
+        if ((w*w-1)*(w*w-1) < EPS(value_t))
+        	return quat<value_t>(0);
+
+        const value_t inv = 1.0/std::sqrt(x*x+y*y+z*z);
+        const value_t fac = std::acos(w)*inv;
 
         return quat<value_t> (0, x*fac, y*fac, z*fac);
     }
@@ -292,7 +295,7 @@ struct dualquat {
         if (x*x+y*y+z*z < EPS(value_t))
         	return dualquat<value_t>(1, 0, 0, 0, 0, X, Y, Z);
 
-        const value_t theta = 2.0*std::sqrt(x*x+y*y+z*z+EPS(value_t));
+        const value_t theta = 2.0*std::sqrt(x*x+y*y+z*z);
         const value_t invvv = 2.0/theta;
         const value_t lx = x*invvv;
         const value_t ly = y*invvv;
@@ -338,11 +341,11 @@ struct dualquat {
 
         assert(isunit());
 
-        if ((w-1)*(w-1) < EPS(value_t))
+        if ((w*w-1)*(w*w-1) < EPS(value_t))
         	return dualquat<value_t>(0, 0, 0, 0, 0, X, Y, Z);
 
-        const value_t theta =  2.0*std::acos(w > 1 ? 1 : w < -1 ? -1 : w);
-        const value_t invvv =  0.5/std::sqrt(x*x+y*y+z*z+EPS(value_t));
+        const value_t theta =  2.0*std::acos(w);
+        const value_t invvv =  0.5/std::sqrt(x*x+y*y+z*z);
         const value_t pitch = -4.0*W*invvv;
         const value_t alpha =  pitch*w;
 

diff --git a/main.cpp b/main.cpp
@@ -2,7 +2,7 @@
 #include <random>
 
 #define NUM_TESTS (1<<22)
-#define THRESHOLD (1E-8)
+#define THRESHOLD (1E-10)
 
 int main () {