Fix displayed 3d coordinates showing gibberish (#23485)

scottshambaugh · web-flow · commit fff2a79ce22a · 2023-07-06T09:43:16.000-07:00
* Rework 3d coordinates mouse hover, now displays coordinate of underlying view pane

Get 3d ortho coordinates working

Get non-1 focal lengths working

Docs

* Make 3D coordinates explicitly call out the backing pane they are on

* Save inverse projection matrix for speed

* Don't use deprecated function

---------

Co-authored-by: Scott Shambaugh &lt;scottshambaugh@users.noreply.github.com&gt;
diff --git a/doc/users/next_whats_new/3d_hover_coordinates.rst b/doc/users/next_whats_new/3d_hover_coordinates.rst
@@ -0,0 +1,10 @@
+3D hover coordinates
+--------------------
+
+The x, y, z coordinates displayed in 3D plots were previously showing
+nonsensical values. This has been fixed to report the coordinate on the view
+pane directly beneath the mouse cursor. This is likely to be most useful when
+viewing 3D plots along a primary axis direction when using an orthographic
+projection, or when a 2D plot has been projected onto one of the 3D axis panes.
+Note that there is still no way to directly display the coordinates of plotted
+data points.
diff --git a/lib/mpl_toolkits/mplot3d/axes3d.py b/lib/mpl_toolkits/mplot3d/axes3d.py
@@ -164,6 +164,7 @@ def __init__(
         # Enable drawing of axes by Axes3D class
         self.set_axis_on()
         self.M = None
+        self.invM = None
 
         # func used to format z -- fall back on major formatters
         self.fmt_zdata = None
@@ -455,6 +456,7 @@ def draw(self, renderer):
 
         # add the projection matrix to the renderer
         self.M = self.get_proj()
+        self.invM = np.linalg.inv(self.M)
 
         collections_and_patches = (
             artist for artist in self._children
@@ -1060,45 +1062,97 @@ def format_zdata(self, z):
             val = func(z)
             return val
 
-    def format_coord(self, xd, yd):
+    def format_coord(self, xv, yv, renderer=None):
         """
-        Given the 2D view coordinates attempt to guess a 3D coordinate.
-        Looks for the nearest edge to the point and then assumes that
-        the point is at the same z location as the nearest point on the edge.
+        Return a string giving the current view rotation angles, or the x, y, z
+        coordinates of the point on the nearest axis pane underneath the mouse
+        cursor, depending on the mouse button pressed.
         """
-
-        if self.M is None:
-            return ''
+        coords = ''
 
         if self.button_pressed in self._rotate_btn:
-            # ignore xd and yd and display angles instead
-            norm_elev = art3d._norm_angle(self.elev)
-            norm_azim = art3d._norm_angle(self.azim)
-            norm_roll = art3d._norm_angle(self.roll)
-            return (f"elevation={norm_elev:.0f}\N{DEGREE SIGN}, "
-                    f"azimuth={norm_azim:.0f}\N{DEGREE SIGN}, "
-                    f"roll={norm_roll:.0f}\N{DEGREE SIGN}"
-                    ).replace("-", "\N{MINUS SIGN}")
-
-        # nearest edge
-        p0, p1 = min(self._tunit_edges(),
-                     key=lambda edge: proj3d._line2d_seg_dist(
-                         (xd, yd), edge[0][:2], edge[1][:2]))
-
-        # scale the z value to match
-        x0, y0, z0 = p0
-        x1, y1, z1 = p1
-        d0 = np.hypot(x0-xd, y0-yd)
-        d1 = np.hypot(x1-xd, y1-yd)
-        dt = d0+d1
-        z = d1/dt * z0 + d0/dt * z1
-
-        x, y, z = proj3d.inv_transform(xd, yd, z, self.M)
-
-        xs = self.format_xdata(x)
-        ys = self.format_ydata(y)
-        zs = self.format_zdata(z)
-        return f'x={xs}, y={ys}, z={zs}'
+            # ignore xv and yv and display angles instead
+            coords = self._rotation_coords()
+
+        elif self.M is not None:
+            coords = self._location_coords(xv, yv, renderer)
+
+        return coords
+
+    def _rotation_coords(self):
+        """
+        Return the rotation angles as a string.
+        """
+        norm_elev = art3d._norm_angle(self.elev)
+        norm_azim = art3d._norm_angle(self.azim)
+        norm_roll = art3d._norm_angle(self.roll)
+        coords = (f"elevation={norm_elev:.0f}\N{DEGREE SIGN}, "
+                  f"azimuth={norm_azim:.0f}\N{DEGREE SIGN}, "
+                  f"roll={norm_roll:.0f}\N{DEGREE SIGN}"
+                  ).replace("-", "\N{MINUS SIGN}")
+        return coords
+
+    def _location_coords(self, xv, yv, renderer):
+        """
+        Return the location on the axis pane underneath the cursor as a string.
+        """
+        p1 = self._calc_coord(xv, yv, renderer)
+        xs = self.format_xdata(p1[0])
+        ys = self.format_ydata(p1[1])
+        zs = self.format_zdata(p1[2])
+        coords = f'x={xs}, y={ys}, z={zs}'
+        return coords
+
+    def _get_camera_loc(self):
+        """
+        Returns the current camera location in data coordinates.
+        """
+        cx, cy, cz, dx, dy, dz = self._get_w_centers_ranges()
+        c = np.array([cx, cy, cz])
+        r = np.array([dx, dy, dz])
+
+        if self._focal_length == np.inf:  # orthographic projection
+            focal_length = 1e9  # large enough to be effectively infinite
+        else:  # perspective projection
+            focal_length = self._focal_length
+        eye = c + self._view_w * self._dist * r / self._box_aspect * focal_length
+        return eye
+
+    def _calc_coord(self, xv, yv, renderer=None):
+        """
+        Given the 2D view coordinates, find the point on the nearest axis pane
+        that lies directly below those coordinates. Returns a 3D point in data
+        coordinates.
+        """
+        if self._focal_length == np.inf:  # orthographic projection
+            zv = 1
+        else:  # perspective projection
+            zv = -1 / self._focal_length
+
+        # Convert point on view plane to data coordinates
+        p1 = np.array(proj3d.inv_transform(xv, yv, zv, self.invM)).ravel()
+
+        # Get the vector from the camera to the point on the view plane
+        vec = self._get_camera_loc() - p1
+
+        # Get the pane locations for each of the axes
+        pane_locs = []
+        for axis in self._axis_map.values():
+            xys, loc = axis.active_pane(renderer)
+            pane_locs.append(loc)
+
+        # Find the distance to the nearest pane by projecting the view vector
+        scales = np.zeros(3)
+        for i in range(3):
+            if vec[i] == 0:
+                scales[i] = np.inf
+            else:
+                scales[i] = (p1[i] - pane_locs[i]) / vec[i]
+        scale = scales[np.argmin(abs(scales))]
+
+        # Calculate the point on the closest pane
+        p2 = p1 - scale*vec
+        return p2
 
     def _on_move(self, event):
         """
@@ -1143,6 +1197,7 @@ def _on_move(self, event):
             self.view_init(elev=elev, azim=azim, roll=roll, share=True)
             self.stale = True
 
+        # Pan
         elif self.button_pressed in self._pan_btn:
             # Start the pan event with pixel coordinates
             px, py = self.transData.transform([self._sx, self._sy])
@@ -1321,21 +1376,27 @@ def _scale_axis_limits(self, scale_x, scale_y, scale_z):
         scale_z : float
             Scale factor for the z data axis.
         """
-        # Get the axis limits and centers
+        # Get the axis centers and ranges
+        cx, cy, cz, dx, dy, dz = self._get_w_centers_ranges()
+
+        # Set the scaled axis limits
+        self.set_xlim3d(cx - dx*scale_x/2, cx + dx*scale_x/2)
+        self.set_ylim3d(cy - dy*scale_y/2, cy + dy*scale_y/2)
+        self.set_zlim3d(cz - dz*scale_z/2, cz + dz*scale_z/2)
+
+    def _get_w_centers_ranges(self):
+        """Get 3D world centers and axis ranges."""
+        # Calculate center of axis limits
         minx, maxx, miny, maxy, minz, maxz = self.get_w_lims()
         cx = (maxx + minx)/2
         cy = (maxy + miny)/2
         cz = (maxz + minz)/2
 
-        # Scale the data range
-        dx = (maxx - minx)*scale_x
-        dy = (maxy - miny)*scale_y
-        dz = (maxz - minz)*scale_z
-
-        # Set the scaled axis limits
-        self.set_xlim3d(cx - dx/2, cx + dx/2)
-        self.set_ylim3d(cy - dy/2, cy + dy/2)
-        self.set_zlim3d(cz - dz/2, cz + dz/2)
+        # Calculate range of axis limits
+        dx = (maxx - minx)
+        dy = (maxy - miny)
+        dz = (maxz - minz)
+        return cx, cy, cz, dx, dy, dz
 
     def set_zlabel(self, zlabel, fontdict=None, labelpad=None, **kwargs):
         """
diff --git a/lib/mpl_toolkits/mplot3d/axis3d.py b/lib/mpl_toolkits/mplot3d/axis3d.py
@@ -295,6 +295,19 @@ def _get_tickdir(self):
         tickdir = np.roll(info_i, -j)[np.roll(tickdirs_base, j)][i]
         return tickdir
 
+    def active_pane(self, renderer):
+        mins, maxs, centers, deltas, tc, highs = self._get_coord_info(renderer)
+        info = self._axinfo
+        index = info['i']
+        if not highs[index]:
+            loc = mins[index]
+            plane = self._PLANES[2 * index]
+        else:
+            loc = maxs[index]
+            plane = self._PLANES[2 * index + 1]
+        xys = np.array([tc[p] for p in plane])
+        return xys, loc
+
     def draw_pane(self, renderer):
         """
         Draw pane.
@@ -304,20 +317,9 @@ def draw_pane(self, renderer):
         renderer : `~matplotlib.backend_bases.RendererBase` subclass
         """
         renderer.open_group('pane3d', gid=self.get_gid())
-
-        mins, maxs, centers, deltas, tc, highs = self._get_coord_info(renderer)
-
-        info = self._axinfo
-        index = info['i']
-        if not highs[index]:
-            plane = self._PLANES[2 * index]
-        else:
-            plane = self._PLANES[2 * index + 1]
-        xys = np.asarray([tc[p] for p in plane])
-        xys = xys[:, :2]
-        self.pane.xy = xys
+        xys, loc = self.active_pane(renderer)
+        self.pane.xy = xys[:, :2]
         self.pane.draw(renderer)
-
         renderer.close_group('pane3d')
 
     @artist.allow_rasterization
diff --git a/lib/mpl_toolkits/mplot3d/proj3d.py b/lib/mpl_toolkits/mplot3d/proj3d.py
@@ -3,34 +3,10 @@
 """
 
 import numpy as np
-import numpy.linalg as linalg
 
 from matplotlib import _api
 
 
-def _line2d_seg_dist(p, s0, s1):
-    """
-    Return the distance(s) from point(s) *p* to segment(s) (*s0*, *s1*).
-
-    Parameters
-    ----------
-    p : (ndim,) or (N, ndim) array-like
-        The points from which the distances are computed.
-    s0, s1 : (ndim,) or (N, ndim) array-like
-        The xy(z...) coordinates of the segment endpoints.
-    """
-    s0 = np.asarray(s0)
-    s01 = s1 - s0  # shape (ndim,) or (N, ndim)
-    s0p = p - s0  # shape (ndim,) or (N, ndim)
-    l2 = s01 @ s01  # squared segment length
-    # Avoid div. by zero for degenerate segments (for them, s01 = (0, 0, ...)
-    # so the value of l2 doesn't matter; this just replaces 0/0 by 0/1).
-    l2 = np.where(l2, l2, 1)
-    # Project onto segment, without going past segment ends.
-    p1 = s0 + np.multiply.outer(np.clip(s0p @ s01 / l2, 0, 1), s01)
-    return ((p - p1) ** 2).sum(axis=-1) ** (1/2)
-
-
 def world_transformation(xmin, xmax,
                          ymin, ymax,
                          zmin, zmax, pb_aspect=None):
@@ -213,17 +189,17 @@ def _proj_transform_vec_clip(vec, M):
     return txs, tys, tzs, tis
 
 
-def inv_transform(xs, ys, zs, M):
+def inv_transform(xs, ys, zs, invM):
     """
-    Transform the points by the inverse of the projection matrix *M*.
+    Transform the points by the inverse of the projection matrix, *invM*.
     """
-    iM = linalg.inv(M)
     vec = _vec_pad_ones(xs, ys, zs)
-    vecr = np.dot(iM, vec)
-    try:
-        vecr = vecr / vecr[3]
-    except OverflowError:
-        pass
+    vecr = np.dot(invM, vec)
+    if vecr.shape == (4,):
+        vecr = vecr.reshape((4, 1))
+    for i in range(vecr.shape[1]):
+        if vecr[3][i] != 0:
+            vecr[:, i] = vecr[:, i] / vecr[3][i]
     return vecr[0], vecr[1], vecr[2]
 
 
diff --git a/lib/mpl_toolkits/mplot3d/tests/baseline_images/test_axes3d/proj3d_lines_dists.png b/lib/mpl_toolkits/mplot3d/tests/baseline_images/test_axes3d/proj3d_lines_dists.png
diff --git a/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py b/lib/mpl_toolkits/mplot3d/tests/test_axes3d.py
@@ -1061,13 +1061,14 @@ def _test_proj_make_M():
 
 def test_proj_transform():
     M = _test_proj_make_M()
+    invM = np.linalg.inv(M)
 
     xs = np.array([0, 1, 1, 0, 0, 0, 1, 1, 0, 0]) * 300.0
     ys = np.array([0, 0, 1, 1, 0, 0, 0, 1, 1, 0]) * 300.0
     zs = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1]) * 300.0
 
     txs, tys, tzs = proj3d.proj_transform(xs, ys, zs, M)
-    ixs, iys, izs = proj3d.inv_transform(txs, tys, tzs, M)
+    ixs, iys, izs = proj3d.inv_transform(txs, tys, tzs, invM)
 
     np.testing.assert_almost_equal(ixs, xs)
     np.testing.assert_almost_equal(iys, ys)
@@ -1154,39 +1155,6 @@ def test_world():
                                 [0, 0, 0, 1]])
 
 
-@mpl3d_image_comparison(['proj3d_lines_dists.png'], style='mpl20')
-def test_lines_dists():
-    fig, ax = plt.subplots(figsize=(4, 6), subplot_kw=dict(aspect='equal'))
-
-    xs = (0, 30)
-    ys = (20, 150)
-    ax.plot(xs, ys)
-    p0, p1 = zip(xs, ys)
-
-    xs = (0, 0, 20, 30)
-    ys = (100, 150, 30, 200)
-    ax.scatter(xs, ys)
-
-    dist0 = proj3d._line2d_seg_dist((xs[0], ys[0]), p0, p1)
-    dist = proj3d._line2d_seg_dist(np.array((xs, ys)).T, p0, p1)
-    assert dist0 == dist[0]
-
-    for x, y, d in zip(xs, ys, dist):
-        c = Circle((x, y), d, fill=0)
-        ax.add_patch(c)
-
-    ax.set_xlim(-50, 150)
-    ax.set_ylim(0, 300)
-
-
-def test_lines_dists_nowarning():
-    # No RuntimeWarning must be emitted for degenerate segments, see GH#22624.
-    s0 = (10, 30, 50)
-    p = (20, 150, 180)
-    proj3d._line2d_seg_dist(p, s0, s0)
-    proj3d._line2d_seg_dist(np.array(p), s0, s0)
-
-
 def test_autoscale():
     fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
     assert ax.get_zscale() == 'linear'
@@ -1963,16 +1931,30 @@ def test_format_coord():
     ax = fig.add_subplot(projection='3d')
     x = np.arange(10)
     ax.plot(x, np.sin(x))
+    xv = 0.1
+    yv = 0.1
     fig.canvas.draw()
-    assert ax.format_coord(0, 0) == 'x=1.8066, y=1.0367, z=−0.0553'
+    assert ax.format_coord(xv, yv) == 'x=10.5227, y=1.0417, z=0.1444'
+
     # Modify parameters
     ax.view_init(roll=30, vertical_axis="y")
     fig.canvas.draw()
-    assert ax.format_coord(0, 0) == 'x=9.1651, y=−0.9215, z=−0.0359'
+    assert ax.format_coord(xv, yv) == 'x=9.1875, y=0.9761, z=0.1291'
+
     # Reset parameters
     ax.view_init()
     fig.canvas.draw()
-    assert ax.format_coord(0, 0) == 'x=1.8066, y=1.0367, z=−0.0553'
+    assert ax.format_coord(xv, yv) == 'x=10.5227, y=1.0417, z=0.1444'
+
+    # Check orthographic projection
+    ax.set_proj_type('ortho')
+    fig.canvas.draw()
+    assert ax.format_coord(xv, yv) == 'x=10.8869, y=1.0417, z=0.1528'
+
+    # Check non-default perspective projection
+    ax.set_proj_type('persp', focal_length=0.1)
+    fig.canvas.draw()
+    assert ax.format_coord(xv, yv) == 'x=9.0620, y=1.0417, z=0.1110'
 
 
 def test_get_axis_position():
