Uniformly distributed random rotations, unit vectors

CHANGELOG.md

@@ -5,8 +5,11 @@ documented here.
 This project adheres to [Semantic Versioning](http://semver.org/).
 
 ## [0.16.0] - WIP
+## Modified
+  * Adjust `UnitQuaternion`s and `Rotation3`s generated from the `Standard` distribution to be uniformly distributed.
 ### Added
   * Add construction of a `Point` from an array by implementing the `From` trait.
+  * Add support for generating uniformly distributed random unit column vectors using the `Standard` distribution.
 
 ## [0.15.0]
 The most notable change of this release is the support for using part of the library without the rust standard

src/base/construction.rs

@@ -6,12 +6,12 @@ use quickcheck::{Arbitrary, Gen};
 use num::{Bounded, One, Zero};
 #[cfg(feature = "std")]
 use rand;
-use rand::distributions::{Distribution, Standard};
+use rand::distributions::{Distribution, Standard, StandardNormal};
 use rand::Rng;
 use std::iter;
 use typenum::{self, Cmp, Greater};
 
-use alga::general::{ClosedAdd, ClosedMul};
+use alga::general::{ClosedAdd, ClosedMul, Real};
 
 use base::allocator::Allocator;
 use base::dimension::{Dim, DimName, Dynamic, U1, U2, U3, U4, U5, U6};
@@ -501,6 +501,19 @@ where
     }
 }
 
+#[cfg(feature = "std")]
+impl<N: Real, D: DimName> Distribution<Unit<VectorN<N, D>>> for Standard
+where
+    DefaultAllocator: Allocator<N, D>,
+    StandardNormal: Distribution<N>,
+{
+    #[inline]
+    fn sample<'a, G: Rng + ?Sized>(&self, rng: &'a mut G) -> Unit<VectorN<N, D>> {
+        Unit::new_normalize(VectorN::from_distribution_generic(D::name(), U1, &mut StandardNormal, rng))
+    }
+}
+
+
 /*
  *
  * Constructors for small matrices and vectors.

src/geometry/quaternion_construction.rs

@@ -6,7 +6,7 @@ use base::storage::Owned;
 use quickcheck::{Arbitrary, Gen};
 
 use num::{One, Zero};
-use rand::distributions::{Distribution, Standard};
+use rand::distributions::{Distribution, Standard, OpenClosed01};
 use rand::Rng;
 
 use alga::general::Real;
@@ -413,11 +413,25 @@ impl<N: Real> One for UnitQuaternion<N> {
 
 impl<N: Real> Distribution<UnitQuaternion<N>> for Standard
 where
-    Standard: Distribution<N>,
+    OpenClosed01: Distribution<N>,
 {
     #[inline]
     fn sample<'a, R: Rng + ?Sized>(&self, rng: &'a mut R) -> UnitQuaternion<N> {
-        UnitQuaternion::from_scaled_axis(rng.gen::<Vector3<N>>() * N::two_pi())
+        // Ken Shoemake's Subgroup Algorithm
+        // Uniform random rotations.
+        // In D. Kirk, editor, Graphics Gems III, pages 124-132. Academic, New York, 1992.
+        let x0 = rng.sample(OpenClosed01);
+        let x1 = rng.sample(OpenClosed01);
+        let x2 = rng.sample(OpenClosed01);
+        let theta1 = N::two_pi() * x1;
+        let theta2 = N::two_pi() * x2;
+        let s1 = theta1.sin();
+        let c1 = theta1.cos();
+        let s2 = theta2.sin();
+        let c2 = theta2.cos();
+        let r1 = (N::one() - x0).sqrt();
+        let r2 = x0.sqrt();
+        Unit::new_unchecked(Quaternion::new(s1 * r1, c1 * r1, s2 * r2, c2 * r2))
     }
 }
 
@@ -433,3 +447,18 @@ where
         UnitQuaternion::from_scaled_axis(axisangle)
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use rand::{self, SeedableRng};
+
+    #[test]
+    fn random_unit_quats_are_unit() {
+        let mut rng = rand::prng::XorShiftRng::from_seed([0xAB; 16]);
+        for _ in 0..1000 {
+            let x = rng.gen::<UnitQuaternion<f32>>();
+            assert!(relative_eq!(x.unwrap().norm(), 1.0))
+        }
+    }
+}

src/geometry/rotation_specialization.rs

@@ -5,7 +5,7 @@ use quickcheck::{Arbitrary, Gen};
 
 use alga::general::Real;
 use num::Zero;
-use rand::distributions::{Distribution, Standard};
+use rand::distributions::{Distribution, Standard, OpenClosed01};
 use rand::Rng;
 use std::ops::Neg;
 
@@ -384,11 +384,34 @@ impl<N: Real> Rotation3<N> {
 
 impl<N: Real> Distribution<Rotation3<N>> for Standard
 where
-    Standard: Distribution<N>,
+    OpenClosed01: Distribution<N>,
 {
     #[inline]
     fn sample<'a, R: Rng + ?Sized>(&self, rng: &mut R) -> Rotation3<N> {
-        Rotation3::new(rng.gen::<Vector3<N>>() * N::two_pi())
+        // James Arvo.
+        // Fast random rotation matrices.
+        // In D. Kirk, editor, Graphics Gems III, pages 117-120. Academic, New York, 1992.
+
+        // Compute a random rotation around Z
+        let theta = N::two_pi() * rng.sample(OpenClosed01);
+        let (ts, tc) = theta.sin_cos();
+        let a = MatrixN::<N, U3>::new(
+            tc, ts, N::zero(),
+            -ts, tc, N::zero(),
+            N::zero(), N::zero(), N::one()
+        );
+
+        // Compute a random rotation *of* Z
+        let phi = N::two_pi() * rng.sample(OpenClosed01);
+        let z = rng.sample(OpenClosed01);
+        let (ps, pc) = phi.sin_cos();
+        let sqrt_z = z.sqrt();
+        let v = Vector3::new(pc * sqrt_z, ps * sqrt_z, (N::one() - z).sqrt());
+        let mut b = v * v.transpose();
+        b += b;
+        b -= MatrixN::<N, U3>::identity();
+
+        Rotation3::from_matrix_unchecked(b * a)
     }
 }