From 352e71656ddb2f49117018f1e0d285fcca7f0372 Mon Sep 17 00:00:00 2001 From: Benjamin Saunders Date: Tue, 3 Jul 2018 19:06:03 -0700 Subject: [PATCH] Uniformly distributed random rotations, unit vectors --- CHANGELOG.md | 3 +++ src/base/construction.rs | 17 ++++++++++-- src/geometry/quaternion_construction.rs | 35 ++++++++++++++++++++++--- src/geometry/rotation_specialization.rs | 29 +++++++++++++++++--- 4 files changed, 76 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0fd586f4..d4ee221b 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/base/construction.rs b/src/base/construction.rs index b6c36735..a78f6a3b 100644 --- a/src/base/construction.rs +++ b/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 Distribution>> for Standard +where + DefaultAllocator: Allocator, + StandardNormal: Distribution, +{ + #[inline] + fn sample<'a, G: Rng + ?Sized>(&self, rng: &'a mut G) -> Unit> { + Unit::new_normalize(VectorN::from_distribution_generic(D::name(), U1, &mut StandardNormal, rng)) + } +} + + /* * * Constructors for small matrices and vectors. diff --git a/src/geometry/quaternion_construction.rs b/src/geometry/quaternion_construction.rs index e18c9948..aa156eef 100644 --- a/src/geometry/quaternion_construction.rs +++ b/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 One for UnitQuaternion { impl Distribution> for Standard where - Standard: Distribution, + OpenClosed01: Distribution, { #[inline] fn sample<'a, R: Rng + ?Sized>(&self, rng: &'a mut R) -> UnitQuaternion { - UnitQuaternion::from_scaled_axis(rng.gen::>() * 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::>(); + assert!(relative_eq!(x.unwrap().norm(), 1.0)) + } + } +} diff --git a/src/geometry/rotation_specialization.rs b/src/geometry/rotation_specialization.rs index 337883a3..e5cd4b4a 100644 --- a/src/geometry/rotation_specialization.rs +++ b/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 Rotation3 { impl Distribution> for Standard where - Standard: Distribution, + OpenClosed01: Distribution, { #[inline] fn sample<'a, R: Rng + ?Sized>(&self, rng: &mut R) -> Rotation3 { - Rotation3::new(rng.gen::>() * 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::::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::::identity(); + + Rotation3::from_matrix_unchecked(b * a) } }