Repository: bitshifter/glam-rs
Pull request #387 (2023-03-16 14:23:08)

Calculate Quat::to_axis_angle in a more numerically stable way
The current implementation results in unnormalized axis when the angle is close to 0.

This way it is done in this PR is probably a bit slower.
But handles angles close to 0 better.
It seems to be how Eigen and nalgebra does it for example.
It guarantees that the axis is normalized.

bitshifter: Not sure why the codegen lint failed since it looks like you have run codegen. Possibly you need to run it again or maybe the version of Rust you are on is producing different codegen.

Otherwise this looks fine in principle.
jnises: > Not sure why the codegen lint failed since it looks like you have run codegen. Possibly you need to run it again or maybe the version of Rust you are on is producing different codegen.
> 
> Otherwise this looks fine in principle.

I hade completely missed the template bit  :)
I was thinking about how error prone it was to update all those files manually 😅 
I still get codegen mismatches when I run that locally, but that is due to line endings it seems. Hopefully they match on CI.

Patch:
diff --git a/codegen/templates/quat.rs.tera b/codegen/templates/quat.rs.tera
index b255396..1686c04 100644
--- a/codegen/templates/quat.rs.tera
+++ b/codegen/templates/quat.rs.tera
@@ -456,18 +456,18 @@ impl {{ self_t }} {
         }
     }
 
-    /// Returns the rotation axis and angle (in radians) of `self`.
+    /// Returns the rotation axis (normalized) and angle (in radians) of `self`.
     #[inline]
     pub fn to_axis_angle(self) -> ({{ vec3_t }}, {{ scalar_t }}) {
         const EPSILON: {{ scalar_t }} = 1.0e-8;
-        const EPSILON_SQUARED: {{ scalar_t }} = EPSILON * EPSILON;
-        let w = self.w;
-        let angle = w.acos_approx() * 2.0;
-        let scale_sq = {{ scalar_t }}::max(1.0 - w * w, 0.0);
-        if scale_sq >= EPSILON_SQUARED {
-            ({{ vec3_t }}::new(self.x, self.y, self.z) * scale_sq.sqrt().recip(), angle)
+        let v = {{ vec3_t }}::new(self.x, self.y, self.z);
+        let length = v.length();
+        if length >= EPSILON {
+            let angle = 2.0 * {{ scalar_t }}::atan2(length, self.w);
+            let axis = v / length;
+            (axis, angle)
         } else {
-            ({{ vec3_t }}::X, angle)
+            ({{ vec3_t }}::X, 0.0)
         }
     }
 
diff --git a/src/f32/coresimd/quat.rs b/src/f32/coresimd/quat.rs
index 6d80256..65454fd 100644
--- a/src/f32/coresimd/quat.rs
+++ b/src/f32/coresimd/quat.rs
@@ -343,21 +343,18 @@ impl Quat {
         }
     }
 
-    /// Returns the rotation axis and angle (in radians) of `self`.
+    /// Returns the rotation axis (normalized) and angle (in radians) of `self`.
     #[inline]
     pub fn to_axis_angle(self) -> (Vec3, f32) {
         const EPSILON: f32 = 1.0e-8;
-        const EPSILON_SQUARED: f32 = EPSILON * EPSILON;
-        let w = self.w;
-        let angle = w.acos_approx() * 2.0;
-        let scale_sq = f32::max(1.0 - w * w, 0.0);
-        if scale_sq >= EPSILON_SQUARED {
-            (
-                Vec3::new(self.x, self.y, self.z) * scale_sq.sqrt().recip(),
-                angle,
-            )
+        let v = Vec3::new(self.x, self.y, self.z);
+        let length = v.length();
+        if length >= EPSILON {
+            let angle = 2.0 * f32::atan2(length, self.w);
+            let axis = v / length;
+            (axis, angle)
         } else {
-            (Vec3::X, angle)
+            (Vec3::X, 0.0)
         }
     }
 
diff --git a/src/f32/scalar/quat.rs b/src/f32/scalar/quat.rs
index 0b72e8f..3ec4f96 100644
--- a/src/f32/scalar/quat.rs
+++ b/src/f32/scalar/quat.rs
@@ -351,21 +351,18 @@ impl Quat {
         }
     }
 
-    /// Returns the rotation axis and angle (in radians) of `self`.
+    /// Returns the rotation axis (normalized) and angle (in radians) of `self`.
     #[inline]
     pub fn to_axis_angle(self) -> (Vec3, f32) {
         const EPSILON: f32 = 1.0e-8;
-        const EPSILON_SQUARED: f32 = EPSILON * EPSILON;
-        let w = self.w;
-        let angle = w.acos_approx() * 2.0;
-        let scale_sq = f32::max(1.0 - w * w, 0.0);
-        if scale_sq >= EPSILON_SQUARED {
-            (
-                Vec3::new(self.x, self.y, self.z) * scale_sq.sqrt().recip(),
-                angle,
-            )
+        let v = Vec3::new(self.x, self.y, self.z);
+        let length = v.length();
+        if length >= EPSILON {
+            let angle = 2.0 * f32::atan2(length, self.w);
+            let axis = v / length;
+            (axis, angle)
         } else {
-            (Vec3::X, angle)
+            (Vec3::X, 0.0)
         }
     }
 
diff --git a/src/f32/sse2/quat.rs b/src/f32/sse2/quat.rs
index 47ad8d7..12fdf84 100644
--- a/src/f32/sse2/quat.rs
+++ b/src/f32/sse2/quat.rs
@@ -350,21 +350,18 @@ impl Quat {
         }
     }
 
-    /// Returns the rotation axis and angle (in radians) of `self`.
+    /// Returns the rotation axis (normalized) and angle (in radians) of `self`.
     #[inline]
     pub fn to_axis_angle(self) -> (Vec3, f32) {
         const EPSILON: f32 = 1.0e-8;
-        const EPSILON_SQUARED: f32 = EPSILON * EPSILON;
-        let w = self.w;
-        let angle = w.acos_approx() * 2.0;
-        let scale_sq = f32::max(1.0 - w * w, 0.0);
-        if scale_sq >= EPSILON_SQUARED {
-            (
-                Vec3::new(self.x, self.y, self.z) * scale_sq.sqrt().recip(),
-                angle,
-            )
+        let v = Vec3::new(self.x, self.y, self.z);
+        let length = v.length();
+        if length >= EPSILON {
+            let angle = 2.0 * f32::atan2(length, self.w);
+            let axis = v / length;
+            (axis, angle)
         } else {
-            (Vec3::X, angle)
+            (Vec3::X, 0.0)
         }
     }
 
diff --git a/src/f32/wasm32/quat.rs b/src/f32/wasm32/quat.rs
index 798faee..e5eed2e 100644
--- a/src/f32/wasm32/quat.rs
+++ b/src/f32/wasm32/quat.rs
@@ -343,21 +343,18 @@ impl Quat {
         }
     }
 
-    /// Returns the rotation axis and angle (in radians) of `self`.
+    /// Returns the rotation axis (normalized) and angle (in radians) of `self`.
     #[inline]
     pub fn to_axis_angle(self) -> (Vec3, f32) {
         const EPSILON: f32 = 1.0e-8;
-        const EPSILON_SQUARED: f32 = EPSILON * EPSILON;
-        let w = self.w;
-        let angle = w.acos_approx() * 2.0;
-        let scale_sq = f32::max(1.0 - w * w, 0.0);
-        if scale_sq >= EPSILON_SQUARED {
-            (
-                Vec3::new(self.x, self.y, self.z) * scale_sq.sqrt().recip(),
-                angle,
-            )
+        let v = Vec3::new(self.x, self.y, self.z);
+        let length = v.length();
+        if length >= EPSILON {
+            let angle = 2.0 * f32::atan2(length, self.w);
+            let axis = v / length;
+            (axis, angle)
         } else {
-            (Vec3::X, angle)
+            (Vec3::X, 0.0)
         }
     }
 
diff --git a/src/f64/dquat.rs b/src/f64/dquat.rs
index 6f89430..66a175c 100644
--- a/src/f64/dquat.rs
+++ b/src/f64/dquat.rs
@@ -341,21 +341,18 @@ impl DQuat {
         }
     }
 
-    /// Returns the rotation axis and angle (in radians) of `self`.
+    /// Returns the rotation axis (normalized) and angle (in radians) of `self`.
     #[inline]
     pub fn to_axis_angle(self) -> (DVec3, f64) {
         const EPSILON: f64 = 1.0e-8;
-        const EPSILON_SQUARED: f64 = EPSILON * EPSILON;
-        let w = self.w;
-        let angle = w.acos_approx() * 2.0;
-        let scale_sq = f64::max(1.0 - w * w, 0.0);
-        if scale_sq >= EPSILON_SQUARED {
-            (
-                DVec3::new(self.x, self.y, self.z) * scale_sq.sqrt().recip(),
-                angle,
-            )
+        let v = DVec3::new(self.x, self.y, self.z);
+        let length = v.length();
+        if length >= EPSILON {
+            let angle = 2.0 * f64::atan2(length, self.w);
+            let axis = v / length;
+            (axis, angle)
         } else {
-            (DVec3::X, angle)
+            (DVec3::X, 0.0)
         }
     }
 


Test patch:
diff --git a/tests/quat.rs b/tests/quat.rs
index 1a6f6e3..499ff40 100644
--- a/tests/quat.rs
+++ b/tests/quat.rs
@@ -507,6 +507,47 @@ macro_rules! impl_quat_tests {
         glam_test!(test_to_array, {
             assert!($new(1.0, 2.0, 3.0, 4.0).to_array() == [1.0, 2.0, 3.0, 4.0]);
         });
+
+        glam_test!(test_to_axis_angle, {
+            {
+                let q = $quat::from_xyzw(
+                    5.28124762e-08,
+                    -5.12559303e-03,
+                    8.29266140e-08,
+                    9.99986828e-01,
+                );
+                assert!(q.is_normalized());
+                let (axis, angle) = q.to_axis_angle();
+                assert!(axis.is_normalized());
+                let q2 = $quat::from_axis_angle(axis, angle);
+                assert!((q.dot(q2) - 1.0).abs() < 1e-6);
+            }
+            {
+                let q = $quat::IDENTITY;
+                let (axis, angle) = q.to_axis_angle();
+                assert!(axis.is_normalized());
+                let q2 = $quat::from_axis_angle(axis, angle);
+                assert!((q.dot(q2) - 1.0).abs() < 1e-6);
+            }
+            {
+                let q = $quat::from_xyzw(0.0, 1.0, 0.0, 0.0);
+                assert!(q.is_normalized());
+                let (axis, angle) = q.to_axis_angle();
+                assert!(axis.is_normalized());
+                let q2 = $quat::from_axis_angle(axis, angle);
+                assert!((q.dot(q2) - 1.0).abs() < 1e-6);
+            }
+            {
+                let axis = $vec3::Z;
+                let angle = core::$t::consts::PI * 0.25;
+                let q = $quat::from_axis_angle(axis, angle);
+                assert!(q.is_normalized());
+                let (axis2, angle2) = q.to_axis_angle();
+                assert!(axis.is_normalized());
+                assert_approx_eq!(axis, axis2);
+                assert_approx_eq!(angle, angle2);
+            }
+        });
     };
 }