diff --git a/manim/utils/space_ops.py b/manim/utils/space_ops.py
index d55c623c5a..450be56c89 100644
--- a/manim/utils/space_ops.py
+++ b/manim/utils/space_ops.py
@@ -151,12 +151,16 @@ def angle_axis_from_quaternion(quaternion: Sequence[float]) -> Sequence[float]:
     Returns
     -------
     Sequence[float]
-        Gives the angle and axis
+        The angle, in the range ``[0, PI]``, and the axis it is measured about.
     """
     axis = normalize(quaternion[1:], fall_back=np.array([1, 0, 0]))
     angle = 2 * np.arccos(quaternion[0])
     if angle > TAU / 2:
+        # A rotation by ``angle`` about ``axis`` is a rotation by
+        # ``TAU - angle`` about ``-axis``, so the axis has to be flipped
+        # along with the angle.
         angle = TAU - angle
+        axis = -axis
     return angle, axis
 
 
diff --git a/tests/module/utils/test_space_ops.py b/tests/module/utils/test_space_ops.py
index bcaf8168ff..09d7af4d4a 100644
--- a/tests/module/utils/test_space_ops.py
+++ b/tests/module/utils/test_space_ops.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+from scipy.spatial.transform import Rotation
 
 from manim.utils.space_ops import *
 from manim.utils.space_ops import shoelace
@@ -54,6 +55,33 @@ def test_rotation_matrices():
     )
 
 
+@pytest.mark.parametrize(
+    "angle",
+    # the first three are at or below PI, where no folding happens
+    [0.5, 2.0, np.pi, 3.5, 4.0, 5.5, 2 * np.pi - 0.01],
+)
+def test_angle_axis_from_quaternion(angle):
+    axis = normalize(np.array([1.0, -2.0, 3.0]))
+    quaternion = quaternion_from_angle_axis(angle, axis)
+    result_angle, result_axis = angle_axis_from_quaternion(quaternion)
+
+    # the returned pair must describe the same rotation as the input
+    np.testing.assert_allclose(
+        rotation_matrix(result_angle, result_axis),
+        rotation_matrix(angle, axis),
+        atol=1e-8,
+    )
+
+    # ... expressed the way scipy's Rotation.as_rotvec expresses it, i.e. with
+    # the angle in [0, PI] and the direction carried by the axis
+    rotation_vector = Rotation.from_quat([*quaternion[1:], quaternion[0]]).as_rotvec()
+    assert 0 <= result_angle <= np.pi
+    np.testing.assert_allclose(result_angle, np.linalg.norm(rotation_vector), atol=1e-8)
+    np.testing.assert_allclose(
+        result_axis, rotation_vector / np.linalg.norm(rotation_vector), atol=1e-8
+    )
+
+
 def test_angle_of_vector():
     assert angle_of_vector(np.array([1, 1, 1])) == np.pi / 4
     assert (