ENH: Add ability to silence Gimbal Lock warnings for rotations.

While it is possible to silence warnings with `catch_warnings()` context
manager and a filter the approach has a few drawbacks:

1) The warnings are still emitted, which can have a perf impact
2) Using catch_warnings hits this bug `python/cpython#73858` that may
   reset a bunch of other warnings.

I also slightly refactor the code by moving the warning 1 level up into
`_compute_euler`
in the stack and having the core function not emit the warnings, but
return a tuple of (data, was_there_a_gimbal_lock). This in particular
also emit the warnings only once per function call instead of once per
gimbal lock.

Another possibility would be to move the warning one more level up into
`as_euler(...)` which is the only public method that uses
`_compute_euler(...)`. The other one being `_as_euler_from_matrix(...)`,
whcih is private and only use in tests. If I can to that, then I think
it might be ok to not add the `gimbal_lock_ok` keyword and just tell
users to reach for the private `_compute_euler`, and discard the warning
flag.

We could also return a boolean array of whcich rotations are Gimbal
Locked. I don't know if that would be useful.

Closes scipy#19306

scipy/spatial/transform/_rotation.pyx

@@ -117,7 +117,7 @@ cdef inline void _quat_canonical(double[:, :] q) noexcept:
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-cdef double[:, :] _compute_euler_from_matrix(
+cdef tuple[double[:, :], bool] _compute_euler_from_matrix(
     np.ndarray[double, ndim=3] matrix, const uchar[:] seq, bint extrinsic=False
 ) noexcept:
     # This is being replaced by the newer: _compute_euler_from_quat
@@ -172,6 +172,8 @@ cdef double[:, :] _compute_euler_from_matrix(
     cdef double eps = 1e-7
     cdef bint safe1, safe2, safe, adjust
 
+    cdef bint gimbal_warning = False
+
     for ind in range(num_rotations):
         _angles = angles[ind, :]
 
@@ -252,15 +254,13 @@ cdef double[:, :] _compute_euler_from_matrix(
 
         # Step 8
         if not safe:
-            warnings.warn("Gimbal lock detected. Setting third angle to zero "
-                          "since it is not possible to uniquely determine "
-                          "all angles.")
+            gimbal_warning = True
 
-    return angles
+    return angles, gimbal_warning
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-cdef double[:, :] _compute_euler_from_quat(
+cdef tuple[double[:, :], bool] _compute_euler_from_quat(
     np.ndarray[double, ndim=2] quat, const uchar[:] seq, bint extrinsic=False
 ) noexcept:
     # The algorithm assumes extrinsic frame transformations. The algorithm
@@ -301,6 +301,8 @@ cdef double[:, :] _compute_euler_from_quat(
     cdef double eps = 1e-7
     cdef int case
 
+    cdef bint gimbal_warning = False
+
     for ind in range(num_rotations):
         _angles = angles[ind, :]
 
@@ -357,11 +359,9 @@ cdef double[:, :] _compute_euler_from_quat(
                 _angles[idx] -= 2 * pi
 
         if case != 0:
-            warnings.warn("Gimbal lock detected. Setting third angle to zero "
-                          "since it is not possible to uniquely determine "
-                          "all angles.")
+            gimbal_warning = True
 
-    return angles
+    return angles, gimbal_warning
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
@@ -1652,7 +1652,7 @@ cdef class Rotation:
             return np.asarray(rotvec)
 
     @cython.embedsignature(True)
-    def _compute_euler(self, seq, degrees, algorithm):
+    def _compute_euler(self, seq, degrees, algorithm, *, gimbal_lock_ok):
         # Prepare axis sequence to call Euler angles conversion algorithm.
 
         if len(seq) != 3:
@@ -1675,25 +1675,28 @@ cdef class Rotation:
             matrix = self.as_matrix()
             if matrix.ndim == 2:
                 matrix = matrix[None, :, :]
-            angles = np.asarray(_compute_euler_from_matrix(
-                matrix, seq.encode(), extrinsic))
+            data, gimbal_warning = _compute_euler_from_matrix(matrix, seq.encode(), extrinsic)
         elif algorithm == 'from_quat':
             quat = self.as_quat()
             if quat.ndim == 1:
                 quat = quat[None, :]
-            angles = np.asarray(_compute_euler_from_quat(
-                    quat, seq.encode(), extrinsic))
+            data, gimbal_warning = _compute_euler_from_quat(quat, seq.encode(), extrinsic)
         else:
             # algorithm can only be 'from_quat' or 'from_matrix'
             assert False
-
+        angles = np.asarray(data)
+        if gimbal_warning and not gimbal_lock_ok:
+            warnings.warn("Gimbal lock detected. Setting third angle to zero "
+                          "since it is not possible to uniquely determine "
+                          "all angles.")
+
         if degrees:
             angles = np.rad2deg(angles)
 
         return angles[0] if self._single else angles
 
     @cython.embedsignature(True)
-    def _as_euler_from_matrix(self, seq, degrees=False):
+    def _as_euler_from_matrix(self, seq, degrees=False,*, gimbal_lock_ok=False):
         """Represent as Euler angles.
 
         Any orientation can be expressed as a composition of 3 elementary
@@ -1720,6 +1723,8 @@ cdef class Rotation:
         degrees : boolean, optional
             Returned angles are in degrees if this flag is True, else they are
             in radians. Default is False.
+        gimbal_lock_ok : boolean, optional
+            Disable the warnings about Gimbal Lock when encountered.
 
         Returns
         -------
@@ -1744,10 +1749,11 @@ cdef class Rotation:
         .. [3] https://en.wikipedia.org/wiki/Gimbal_lock#In_applied_mathematics
 
         """
-        return self._compute_euler(seq, degrees, 'from_matrix')
+        return self._compute_euler(seq, degrees, 'from_matrix',
+                                   gimbal_lock_ok=gimbal_lock_ok)
 
     @cython.embedsignature(True)
-    def as_euler(self, seq, degrees=False):
+    def as_euler(self, seq, degrees=False, *, gimbal_lock_ok=False):
         """Represent as Euler angles.
 
         Any orientation can be expressed as a composition of 3 elementary
@@ -1774,6 +1780,8 @@ cdef class Rotation:
         degrees : boolean, optional
             Returned angles are in degrees if this flag is True, else they are
             in radians. Default is False.
+        gimbal_lock_ok : boolean, optional
+            Disable the warnings about Gimbal Lock when encountered.
 
         Returns
         -------
@@ -1832,8 +1840,17 @@ cdef class Rotation:
         >>> r.as_euler('zxy', degrees=True).shape
         (3, 3)
 
+        Notes
+        -----
+
+        While it is possible to supress the Gimbal Lock warning with
+        ``catch_warnings`` and ``filters``, this still has performance impact
+        and create other issues with problems:
+
+        - https://github.com/python/cpython/issues/73858
+
         """
-        return self._compute_euler(seq, degrees, 'from_quat')
+        return self._compute_euler(seq, degrees, 'from_quat', gimbal_lock_ok=gimbal_lock_ok)
 
     @cython.embedsignature(True)
     def as_mrp(self):

scipy/spatial/transform/tests/test_rotation.py

@@ -9,6 +9,8 @@
 
 import pickle
 import copy
+from contextlib import contextmanager
+import warnings
 
 def test_generic_quat_matrix():
     x = np.array([[3, 4, 0, 0], [5, 12, 0, 0]])
@@ -683,9 +685,22 @@ def test_stats(error, mean_max, rms_max):
     test_stats(angles_mat - angles, 1e-15, 1e-13)
 
 
+@contextmanager
+def maybe_warn_gimbal_lock(should_warn):
+    if should_warn:
+        with pytest.warns(UserWarning, match="Gimbal lock"):
+            yield
+
+    else:
+        with warnings.catch_warnings():
+            warnings.simplefilter("error")
+            yield
+
+
 @pytest.mark.parametrize("seq_tuple", permutations("xyz"))
 @pytest.mark.parametrize("intrinsic", (False, True))
-def test_as_euler_degenerate_asymmetric_axes(seq_tuple, intrinsic):
+@pytest.mark.parametrize("gimbal_lock_ok", (False, True))
+def test_as_euler_degenerate_asymmetric_axes(seq_tuple, intrinsic, gimbal_lock_ok):
     # Since we cannot check for angle equality, we check for rotation matrix
     # equality
     angles = np.array([
@@ -700,16 +715,19 @@ def test_as_euler_degenerate_asymmetric_axes(seq_tuple, intrinsic):
     rotation = Rotation.from_euler(seq, angles, degrees=True)
     mat_expected = rotation.as_matrix()
 
-    with pytest.warns(UserWarning, match="Gimbal lock"):
-        angle_estimates = rotation.as_euler(seq, degrees=True)
+    with maybe_warn_gimbal_lock(not gimbal_lock_ok):
+        angle_estimates = rotation.as_euler(
+            seq, degrees=True, gimbal_lock_ok=gimbal_lock_ok
+        )
     mat_estimated = Rotation.from_euler(seq, angle_estimates, degrees=True).as_matrix()
 
     assert_array_almost_equal(mat_expected, mat_estimated)
 
 
 @pytest.mark.parametrize("seq_tuple", permutations("xyz"))
 @pytest.mark.parametrize("intrinsic", (False, True))
-def test_as_euler_degenerate_symmetric_axes(seq_tuple, intrinsic):
+@pytest.mark.parametrize("gimbal_lock_ok", (False, True))
+def test_as_euler_degenerate_symmetric_axes(seq_tuple, intrinsic, gimbal_lock_ok):
     # Since we cannot check for angle equality, we check for rotation matrix
     # equality
     angles = np.array([
@@ -724,16 +742,19 @@ def test_as_euler_degenerate_symmetric_axes(seq_tuple, intrinsic):
     rotation = Rotation.from_euler(seq, angles, degrees=True)
     mat_expected = rotation.as_matrix()
 
-    with pytest.warns(UserWarning, match="Gimbal lock"):
-        angle_estimates = rotation.as_euler(seq, degrees=True)
+    with maybe_warn_gimbal_lock(not gimbal_lock_ok):
+        angle_estimates = rotation.as_euler(
+            seq, degrees=True, gimbal_lock_ok=gimbal_lock_ok
+        )
     mat_estimated = Rotation.from_euler(seq, angle_estimates, degrees=True).as_matrix()
 
     assert_array_almost_equal(mat_expected, mat_estimated)
 
 
 @pytest.mark.parametrize("seq_tuple", permutations("xyz"))
 @pytest.mark.parametrize("intrinsic", (False, True))
-def test_as_euler_degenerate_compare_algorithms(seq_tuple, intrinsic):
+@pytest.mark.parametrize("gimbal_lock_ok", (False, True))
+def test_as_euler_degenerate_compare_algorithms(seq_tuple, intrinsic, gimbal_lock_ok):
     # this test makes sure that both algorithms are doing the same choices
     # in degenerate cases
 
@@ -751,8 +772,8 @@ def test_as_euler_degenerate_compare_algorithms(seq_tuple, intrinsic):
     rot = Rotation.from_euler(seq, angles, degrees=True)
     with pytest.warns(UserWarning, match="Gimbal lock"):
         estimates_matrix = rot._as_euler_from_matrix(seq, degrees=True)
-    with pytest.warns(UserWarning, match="Gimbal lock"):
-        estimates_quat = rot.as_euler(seq, degrees=True)
+    with maybe_warn_gimbal_lock(not gimbal_lock_ok):
+        estimates_quat = rot.as_euler(seq, degrees=True, gimbal_lock_ok=gimbal_lock_ok)
     assert_allclose(
         estimates_matrix[:, [0, 2]], estimates_quat[:, [0, 2]], atol=0, rtol=1e-12
     )
@@ -778,8 +799,8 @@ def test_as_euler_degenerate_compare_algorithms(seq_tuple, intrinsic):
     rot = Rotation.from_euler(seq, angles, degrees=True)
     with pytest.warns(UserWarning, match="Gimbal lock"):
         estimates_matrix = rot._as_euler_from_matrix(seq, degrees=True)
-    with pytest.warns(UserWarning, match="Gimbal lock"):
-        estimates_quat = rot.as_euler(seq, degrees=True)
+    with maybe_warn_gimbal_lock(not gimbal_lock_ok):
+        estimates_quat = rot.as_euler(seq, degrees=True, gimbal_lock_ok=gimbal_lock_ok)
     assert_allclose(
         estimates_matrix[:, [0, 2]], estimates_quat[:, [0, 2]], atol=0, rtol=1e-12
     )