Natural 3D rotation with mouse

- Addresses Issue matplotlib#28288
- Introduces three-dimensional rotation by mouse using a variation on Ken Shoemake's ARCBALL
- Provides a minimal Quaternion class, to avoid an additional  dependency on a large package like 'numpy-quaternion'

.gitignore

@@ -114,3 +114,4 @@ lib/matplotlib/backends/web_backend/node_modules/
 lib/matplotlib/backends/web_backend/package-lock.json
 
 LICENSE/LICENSE_QHULL
+/.venv

doc/users/next_whats_new/mouse_rotation.rst

@@ -0,0 +1,12 @@
+Rotating 3d plots with the mouse
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Rotating three-dimensional plots with the mouse has been made more intuitive.
+The plot now reacts the same way to mouse movement, independent of the
+particular orientation at hand; and it is possible to control all 3 rotational
+degrees of freedom (azimuth, elevation, and roll). It uses a variation on
+Ken Shoemake's ARCBALL [Shoemake1992]_.
+
+.. [Shoemake1992] Ken Shoemake, "ARCBALL: A user interface for specifying
+  three-dimensional rotation using a mouse." in Proceedings of Graphics
+  Interface '92, 1992, pp. 151-156, https://doi.org/10.20380/GI1992.18

lib/mpl_toolkits/mplot3d/axes3d.py

@@ -1484,6 +1484,19 @@ def _calc_coord(self, xv, yv, renderer=None):
         p2 = p1 - scale*vec
         return p2, pane_idx
 
+    def _arcball(self, p):
+        """
+        Convert a point p = [x, y] to a point on the virtual trackball
+        """
+        p = 2 * p
+        r = p[0]**2 + p[1]**2
+        # Ken Shoemake's arcball
+        if r > 1:
+            p = np.concatenate(([0], p/math.sqrt(r)))
+        else:
+            p = np.concatenate(([math.sqrt(1-r)], p))
+        return p
+
     def _on_move(self, event):
         """
         Mouse moving.
@@ -1519,11 +1532,26 @@ def _on_move(self, event):
             if dx == 0 and dy == 0:
                 return
 
+            # Convert to quaternion
+            elev = np.deg2rad(self.elev)
+            azim = np.deg2rad(self.azim)
             roll = np.deg2rad(self.roll)
-            delev = -(dy/h)*180*np.cos(roll) + (dx/w)*180*np.sin(roll)
-            dazim = -(dy/h)*180*np.sin(roll) - (dx/w)*180*np.cos(roll)
-            elev = self.elev + delev
-            azim = self.azim + dazim
+            q = Quaternion.from_cardan_angles(elev, azim, roll)
+
+            # Update quaternion
+            current_point = np.array([self._sx/w, self._sy/h])
+            new_point = np.array([x/w, y/h])
+            current_vec = self._arcball(current_point)
+            new_vec = self._arcball(new_point)
+            dq = Quaternion.from_to(current_vec, new_vec)
+            q = dq*q
+
+            # Convert to elev, azim, roll
+            elev, azim, roll = q.as_cardan_angles()
+            azim = np.rad2deg(azim)
+            elev = np.rad2deg(elev)
+            roll = np.rad2deg(roll)
+
             vertical_axis = self._axis_names[self._vertical_axis]
             self.view_init(
                 elev=elev,
@@ -3706,3 +3734,77 @@ def get_test_data(delta=0.05):
     Y = Y * 10
     Z = Z * 500
     return X, Y, Z
+
+
+class Quaternion:
+    """
+    Quaternions
+    consisting of scalar, along 1, and vector, with components along i, j, k
+    """
+
+    def __init__(self, scalar, vector):
+        self.scalar = scalar
+        self.vector = np.array(vector)
+
+    def __neg__(self):
+        return self.__class__(-self.scalar, -self.vector)
+
+    def __mul__(self, other):
+        """
+        Product of two quaternions
+        i*i = j*j = k*k = i*j*k = -1
+        """
+        return self.__class__(
+            self.scalar*other.scalar - np.dot(self.vector, other.vector),
+            self.scalar*other.vector + self.vector*other.scalar
+            + np.cross(self.vector, other.vector))
+
+    def __eq__(self, other):
+        return (self.scalar == other.scalar) and (self.vector == other.vector).all
+
+    def __repr__(self):
+        return "Quaternion({}, {})".format(repr(self.scalar), repr(self.vector))
+
+    @classmethod
+    def from_to(cls, r1, r2):
+        """The quaternion that rotates vector r1 to vector r2."""
+        k = np.cross(r1, r2)
+        nk = np.linalg.norm(k)
+        th = np.arctan2(nk, np.dot(r1, r2))
+        th = th/2
+        if nk == 0:  # r1 and r2 are parallel or anti-parallel
+            if np.dot(r1, r2) < 0:
+                q = cls(0, [1, 0, 0])  # axis of the circle that r1 and r2 are on
+            else:
+                q = cls(1, [0, 0, 0])  # = 1, no rotation
+        else:
+            q = cls(math.cos(th), k*math.sin(th)/nk)
+        return q
+
+    @classmethod
+    def from_cardan_angles(cls, elev, azim, roll):
+        """
+        Converts the angles to a quaternion
+            q = exp((roll/2)*e_x)*exp((elev/2)*e_y)*exp((-azim/2)*e_z)
+        i.e., the angles are a kind of Tait-Bryan angles, -z,y',x".
+        """
+        ca, sa = np.cos(azim/2), np.sin(azim/2)
+        ce, se = np.cos(elev/2), np.sin(elev/2)
+        cr, sr = np.cos(roll/2), np.sin(roll/2)
+
+        qw = ca*ce*cr + sa*se*sr
+        qx = ca*ce*sr - sa*se*cr
+        qy = ca*se*cr + sa*ce*sr
+        qz = ca*se*sr - sa*ce*cr
+        return cls(qw, [qx, qy, qz])
+
+    def as_cardan_angles(self):
+        """The inverse of from_cardan_angles()."""
+        qw = self.scalar
+        qx = self.vector[..., 0]
+        qy = self.vector[..., 1]
+        qz = self.vector[..., 2]
+        azim = np.arctan2(2*(-qw*qz+qx*qy), qw*qw+qx*qx-qy*qy-qz*qz)
+        elev = np.arcsin( 2*( qw*qy+qz*qx)/(qw*qw+qx*qx+qy*qy+qz*qz))  # noqa E201
+        roll = np.arctan2(2*( qw*qx-qy*qz), qw*qw-qx*qx-qy*qy+qz*qz)   # noqa E201
+        return elev, azim, roll

lib/mpl_toolkits/mplot3d/tests/test_axes3d.py

@@ -1766,6 +1766,82 @@ def test_shared_axes_retick():
     assert ax2.get_zlim() == (-0.5, 2.5)
 
 
+def test_quaternion():
+    # 1:
+    q1 = axes3d.Quaternion(1, [0, 0, 0])
+    assert q1.scalar == 1
+    assert (q1.vector == [0, 0, 0]).all
+    # __neg__:
+    assert (-q1).scalar == -1
+    assert ((-q1).vector == [0, 0, 0]).all
+    # i, j, k:
+    qi = axes3d.Quaternion(0, [1, 0, 0])
+    assert qi.scalar == 0
+    assert (qi.vector == [1, 0, 0]).all
+    qj = axes3d.Quaternion(0, [0, 1, 0])
+    assert qj.scalar == 0
+    assert (qj.vector == [0, 1, 0]).all
+    qk = axes3d.Quaternion(0, [0, 0, 1])
+    assert qk.scalar == 0
+    assert (qk.vector == [0, 0, 1]).all
+    # i^2 = j^2 = k^2 = -1:
+    assert qi*qi == -q1
+    assert qj*qj == -q1
+    assert qk*qk == -q1
+    # identity:
+    assert q1*qi == qi
+    assert q1*qj == qj
+    assert q1*qk == qk
+    # i*j=k, j*k=i, k*i=j:
+    assert qi*qj == qk
+    assert qj*qk == qi
+    assert qk*qi == qj
+    assert qj*qi == -qk
+    assert qk*qj == -qi
+    assert qi*qk == -qj
+    # __mul__:
+    assert (axes3d.Quaternion(2, [3, 4, 5]) * axes3d.Quaternion(6, [7, 8, 9])
+            == axes3d.Quaternion(-86, [28, 48, 44]))
+    # from_to():
+    for r1, r2, q in [
+        ([1, 0, 0], [0, 1, 0], axes3d.Quaternion(np.sqrt(1/2), [0, 0, np.sqrt(1/2)])),
+        ([1, 0, 0], [0, 0, 1], axes3d.Quaternion(np.sqrt(1/2), [0, -np.sqrt(1/2), 0])),
+        ([1, 0, 0], [1, 0, 0], axes3d.Quaternion(1, [0, 0, 0]))
+    ]:
+        assert axes3d.Quaternion.from_to(r1, r2) == q
+    # from_cardan_angles(), as_cardan_angles():
+    for elev, azim, roll in [(0, 0, 0),
+                             (90, 0, 0), (0, 90, 0), (0, 0, 90),
+                             (0, 30, 30), (30, 0, 30), (30, 30, 0)]:
+        q = axes3d.Quaternion.from_cardan_angles(
+            np.deg2rad(elev), np.deg2rad(azim), np.deg2rad(roll))
+        e, a, r = np.rad2deg(axes3d.Quaternion.as_cardan_angles(q))
+        assert np.isclose(e, elev)
+        assert np.isclose(a, azim)
+        assert np.isclose(r, roll)
+
+
+def test_rotate():
+    """Test rotating using the left mouse button."""
+    for roll in [0, 30]:
+        fig = plt.figure()
+        ax = fig.add_subplot(1, 1, 1, projection='3d')
+        ax.view_init(0, 0, roll)
+        ax.figure.canvas.draw()
+
+        # drag mouse horizontally to change orientation
+        ax._button_press(
+            mock_event(ax, button=MouseButton.LEFT, xdata=0, ydata=0))
+        ax._on_move(
+            mock_event(ax, button=MouseButton.LEFT,
+                           xdata=0.5*ax._pseudo_w, ydata=0*ax._pseudo_h))
+        ax.figure.canvas.draw()
+
+        assert np.isclose(ax.elev, roll)
+        assert np.isclose(ax.azim, -90)
+        assert np.isclose(ax.roll, 0)
+
+
 def test_pan():
     """Test mouse panning using the middle mouse button."""