use shoemake algorithm

Former-commit-id: f7766f3b16dfe80aec014125bfa178bfec5524ec

include/igl/random_quaternion.cpp

@@ -10,11 +10,12 @@
 template <typename Scalar>
 IGL_INLINE Eigen::Quaternion<Scalar> igl::random_quaternion()
 {
-  // http://mathproofs.blogspot.com/2005/05/uniformly-distributed-random-unit.html
   const auto & unit_rand = []()->Scalar
   {
     return ((Scalar)rand() / (Scalar)RAND_MAX);
   };
+#ifdef false
+  // http://mathproofs.blogspot.com/2005/05/uniformly-distributed-random-unit.html
   const Scalar t0 = 2.*M_PI*unit_rand();
   const Scalar t1 = acos(1.-2.*unit_rand());
   const Scalar t2 = 0.5*(M_PI*unit_rand() + acos(unit_rand()));
@@ -23,6 +24,57 @@ IGL_INLINE Eigen::Quaternion<Scalar> igl::random_quaternion()
     1.*cos(t0)*sin(t1)*sin(t2),
     1.*cos(t1)*sin(t2),
     1.*cos(t2));
+#elif false
+  // "Uniform Random Rotations" [Shoemake 1992] method 1
+  const auto & uurand = [&unit_rand]()->Scalar
+  {
+    return unit_rand()*2.-1.; 
+  };
+  Scalar x = uurand();
+  Scalar y = uurand();
+  Scalar z = uurand();
+  Scalar w = uurand();
+  const auto & hype = [&uurand](Scalar & x, Scalar & y)->Scalar
+  {
+    Scalar s1;
+    while((s1 = x*x + y*y) > 1.0)
+    {
+      x = uurand();
+      y = uurand();
+    }
+    return s1;
+  };
+  Scalar s1 = hype(x,y);
+  Scalar s2 = hype(z,w);
+  Scalar num1 = -2.*log(s1);
+  Scalar num2 = -2.*log(s2);
+  Scalar r = num1 + num2;
+  Scalar root1 = sqrt((num1/s1)/r);
+  Scalar root2 = sqrt((num2/s2)/r);
+  return Eigen::Quaternion<Scalar>(
+    x*root1,
+    y*root1,
+    z*root2,
+    w*root2);
+#else
+  // Shoemake method 2
+  const Scalar x0 = unit_rand();
+  const Scalar x1 = unit_rand();
+  const Scalar x2 = unit_rand();
+  const Scalar r1 = sqrt(1.0 - x0);
+  const Scalar r2 = sqrt(x0);
+  const Scalar t1 = 2.*M_PI*x1;
+  const Scalar t2 = 2.*M_PI*x2;
+  const Scalar c1 = cos(t1);
+  const Scalar s1 = sin(t1);
+  const Scalar c2 = cos(t2);
+  const Scalar s2 = sin(t2);
+  return Eigen::Quaternion<Scalar>(
+    s1*r1,
+    c1*r1,
+    s2*r2,
+    c2*r2);
+#endif
 }
 
 #ifdef IGL_STATIC_LIBRARY