ENH: allow image to interpolate post RGBA

Co-authored-by: Elliott Sales de Andrade <[email protected]>

Author: jklymak

doc/users/next_whats_new/image_interpolation_rgbastage.rst

@@ -0,0 +1,21 @@
+Image interpolation now possible at RGBA stage
+----------------------------------------------
+
+Images in Matplotlib created via `~.axes.Axes.imshow` are resampled to match 
+the resolution of the current canvas.  It is useful to apply an anto-aliasing
+filter when downsampling to reduce Moire effects.  By default, interpolation 
+is done on the data, a norm applied, and then the colormapping performed. 
+
+However, it is often desireable for the anti-aliasing interpolation to happen 
+in RGBA space, where the colors are interpolated rather than the data.  This 
+usually leads to colors outside the colormap, but visually blends adjacent 
+colors, and is what browsers and other image processing software does. 
+
+A new keyword argument *interpolation_stage* is provided for 
+`~.axes.Axes.imshow` to set the stage at which the anti-aliasing interpolation 
+happens.  The default is the current behaviour of "data", with the alternative
+being "rgba" for the newly-available behavior.  
+
+For more details see the discussion of the new keyword argument in
+:doc:`/gallery/images_contours_and_fields/image_antialiasing`.
+

examples/images_contours_and_fields/image_antialiasing.py

@@ -5,12 +5,21 @@
 
 Images are represented by discrete pixels, either on the screen or in an
 image file.  When data that makes up the image has a different resolution
-than its representation on the screen we will see aliasing effects.
+than its representation on the screen we will see aliasing effects.  How
+noticeable these are depends on how much down-sampling takes place in
+the change of resolution (if any).
 
-The default image interpolation in Matplotlib is 'antialiased'.  This uses a
-hanning interpolation for reduced aliasing in most situations. Only when there
-is upsampling by a factor of 1, 2 or >=3 is 'nearest' neighbor interpolation
-used.
+When subsampling data, aliasing is reduced by smoothing first and then
+subsampling the smoothed data.  In Matplotlib, we can do that
+smoothing before mapping the data to colors, or we can do the smoothing
+on the RGB(A) data in the final image.  The difference between these is
+shown below, and controlled with the *interpolation_stage* keyword argument.
+
+The default image interpolation in Matplotlib is 'antialiased', and
+it is applied to the data.  This uses a
+hanning interpolation on the data provided by the user for reduced aliasing
+in most situations. Only when there is upsampling by a factor of 1, 2 or
+>=3 is 'nearest' neighbor interpolation used.
 
 Other anti-aliasing filters can be specified in `.Axes.imshow` using the
 *interpolation* keyword argument.
@@ -20,26 +29,55 @@
 import matplotlib.pyplot as plt
 
 ###############################################################################
-# First we generate a 500x500 px image with varying frequency content:
-x = np.arange(500) / 500 - 0.5
-y = np.arange(500) / 500 - 0.5
+# First we generate a 450x450 pixel image with varying frequency content:
+N = 450
+x = np.arange(N) / N - 0.5
+y = np.arange(N) / N - 0.5
+aa = np.ones((N, N))
+aa[::2, :] = -1
 
 X, Y = np.meshgrid(x, y)
 R = np.sqrt(X**2 + Y**2)
-f0 = 10
-k = 250
+f0 = 5
+k = 100
 a = np.sin(np.pi * 2 * (f0 * R + k * R**2 / 2))
-
-
+# make the left hand side of this
+a[:int(N / 2), :][R[:int(N / 2), :] < 0.4] = -1
+a[:int(N / 2), :][R[:int(N / 2), :] < 0.3] = 1
+aa[:, int(N / 3):] = a[:, int(N / 3):]
+a = aa
 ###############################################################################
-# The following images are subsampled from 500 data pixels to 303 rendered
-# pixels. The Moire patterns in the 'nearest' interpolation are caused by the
-# high-frequency data being subsampled.  The 'antialiased' image
+# The following images are subsampled from 450 data pixels to either
+# 125 pixels or 250 pixels (depending on your display).
+# The Moire patterns in the 'nearest' interpolation are caused by the
+# high-frequency data being subsampled.  The 'antialiased' imaged
 # still has some Moire patterns as well, but they are greatly reduced.
-fig, axs = plt.subplots(1, 2, figsize=(7, 4), constrained_layout=True)
-for ax, interp in zip(axs, ['nearest', 'antialiased']):
-    ax.imshow(a, interpolation=interp, cmap='gray')
-    ax.set_title(f"interpolation='{interp}'")
+#
+# There are substantial differences between the 'data' interpolation and
+# the 'rgba' interpolation.  The alternating bands of red and blue on the
+# left third of the image are subsampled.  By interpolating in 'data' space
+# (the default) the antialiasing filter makes the stripes close to white,
+# because the average of -1 and +1 is zero, and zero is white in this
+# colormap.
+#
+# Conversely, when the anti-aliasing occurs in 'rgba' space, the red and
+# blue are combined visually to make purple.  This behaviour is more like a
+# typical image processing package, but note that purple is not in the
+# original colormap, so it is no longer possible to invert individual
+# pixels back to their data value.
+
+fig, axs = plt.subplots(2, 2, figsize=(5, 6), constrained_layout=True)
+axs[0, 0].imshow(a, interpolation='nearest', cmap='RdBu_r')
+axs[0, 0].set_xlim(100, 200)
+axs[0, 0].set_ylim(275, 175)
+axs[0, 0].set_title('Zoom')
+
+for ax, interp, space in zip(axs.flat[1:],
+                             ['nearest', 'antialiased', 'antialiased'],
+                             ['data', 'data', 'rgba']):
+    ax.imshow(a, interpolation=interp, interpolation_stage=space,
+              cmap='RdBu_r')
+    ax.set_title(f"interpolation='{interp}'\nspace='{space}'")
 plt.show()
 
 ###############################################################################
@@ -63,7 +101,7 @@
 plt.show()
 
 ###############################################################################
-# Apart from the default 'hanning' antialiasing  `~.Axes.imshow` supports a
+# Apart from the default 'hanning' antialiasing, `~.Axes.imshow` supports a
 # number of different interpolation algorithms, which may work better or
 # worse depending on the pattern.
 fig, axs = plt.subplots(1, 2, figsize=(7, 4), constrained_layout=True)
@@ -72,7 +110,6 @@
     ax.set_title(f"interpolation='{interp}'")
 plt.show()
 
-
 #############################################################################
 #
 # .. admonition:: References

lib/matplotlib/artist.py

@@ -1450,6 +1450,7 @@ def aliased_name(self, s):
         'matplotlib.image._ImageBase.set_filternorm',
         'matplotlib.image._ImageBase.set_filterrad',
         'matplotlib.image._ImageBase.set_interpolation',
+        'matplotlib.image._ImageBase.set_interpolation_stage',
         'matplotlib.image._ImageBase.set_resample',
         'matplotlib.text._AnnotationBase.set_annotation_clip',
     }

lib/matplotlib/axes/_axes.py

@@ -5289,8 +5289,9 @@ def fill_betweenx(self, y, x1, x2=0, where=None,
     @_api.make_keyword_only("3.5", "aspect")
     @_preprocess_data()
     def imshow(self, X, cmap=None, norm=None, aspect=None,
-               interpolation=None, alpha=None, vmin=None, vmax=None,
-               origin=None, extent=None, *, filternorm=True, filterrad=4.0,
+               interpolation=None, alpha=None,
+               vmin=None, vmax=None, origin=None, extent=None, *,
+               interpolation_stage=None, filternorm=True, filterrad=4.0,
                resample=None, url=None, **kwargs):
         """
         Display data as an image, i.e., on a 2D regular raster.
@@ -5382,6 +5383,12 @@ def imshow(self, X, cmap=None, norm=None, aspect=None,
             which can be set by *filterrad*. Additionally, the antigrain image
             resize filter is controlled by the parameter *filternorm*.
 
+        interpolation_stage : {'data', 'rgba'}, default: 'data'
+            If 'data', interpolation
+            is carried out on the data provided by the user.  If 'rgba', the
+            interpolation is carried out after the colormapping has been
+            applied (visual interpolation).
+
         alpha : float or array-like, optional
             The alpha blending value, between 0 (transparent) and 1 (opaque).
             If *alpha* is an array, the alpha blending values are applied pixel
@@ -5482,9 +5489,11 @@ def imshow(self, X, cmap=None, norm=None, aspect=None,
         if aspect is None:
             aspect = rcParams['image.aspect']
         self.set_aspect(aspect)
-        im = mimage.AxesImage(self, cmap, norm, interpolation, origin, extent,
-                              filternorm=filternorm, filterrad=filterrad,
-                              resample=resample, **kwargs)
+        im = mimage.AxesImage(self, cmap, norm, interpolation,
+                              origin, extent, filternorm=filternorm,
+                              filterrad=filterrad, resample=resample,
+                              interpolation_stage=interpolation_stage,
+                              **kwargs)
 
         im.set_data(X)
         im.set_alpha(alpha)

lib/matplotlib/image.py

@@ -238,6 +238,8 @@ def __init__(self, ax,
                  filternorm=True,
                  filterrad=4.0,
                  resample=False,
+                 *,
+                 interpolation_stage=None,
                  **kwargs
                  ):
         martist.Artist.__init__(self)
@@ -249,6 +251,7 @@ def __init__(self, ax,
         self.set_filternorm(filternorm)
         self.set_filterrad(filterrad)
         self.set_interpolation(interpolation)
+        self.set_interpolation_stage(interpolation_stage)
         self.set_resample(resample)
         self.axes = ax
 
@@ -392,8 +395,7 @@ def _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification=1.0,
         if not unsampled:
             if not (A.ndim == 2 or A.ndim == 3 and A.shape[-1] in (3, 4)):
                 raise ValueError(f"Invalid shape {A.shape} for image data")
-
-            if A.ndim == 2:
+            if A.ndim == 2 and self._interpolation_stage != 'rgba':
                 # if we are a 2D array, then we are running through the
                 # norm + colormap transformation.  However, in general the
                 # input data is not going to match the size on the screen so we
@@ -541,6 +543,9 @@ def _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification=1.0,
                                        ):
                     output = self.norm(resampled_masked)
             else:
+                if A.ndim == 2:  # _interpolation_stage == 'rgba'
+                    self.norm.autoscale_None(A)
+                    A = self.to_rgba(A)
                 if A.shape[2] == 3:
                     A = _rgb_to_rgba(A)
                 alpha = self._get_scalar_alpha()
@@ -773,6 +778,22 @@ def set_interpolation(self, s):
         self._interpolation = s
         self.stale = True
 
+    def set_interpolation_stage(self, s):
+        """
+        Set when interpolation happens during the transform to RGBA.
+
+        Parameters
+        ----------
+        s : {'data', 'rgba'} or None
+            Whether to apply up/downsampling interpolation in data or rgba
+            space.
+        """
+        if s is None:
+            s = "data"  # placeholder for maybe having rcParam
+        _api.check_in_list(['data', 'rgba'])
+        self._interpolation_stage = s
+        self.stale = True
+
     def can_composite(self):
         """Return whether the image can be composited with its neighbors."""
         trans = self.get_transform()
@@ -854,6 +875,11 @@ class AxesImage(_ImageBase):
         'bicubic', 'spline16', 'spline36', 'hanning', 'hamming', 'hermite',
         'kaiser', 'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell',
         'sinc', 'lanczos', 'blackman'.
+    interpolation_stage : {'data', 'rgba'}, default: 'data'
+        If 'data', interpolation
+        is carried out on the data provided by the user.  If 'rgba', the
+        interpolation is carried out after the colormapping has been
+        applied (visual interpolation).
     origin : {'upper', 'lower'}, default: :rc:`image.origin`
         Place the [0, 0] index of the array in the upper left or lower left
         corner of the axes. The convention 'upper' is typically used for
@@ -890,6 +916,8 @@ def __init__(self, ax,
                  filternorm=True,
                  filterrad=4.0,
                  resample=False,
+                 *,
+                 interpolation_stage=None,
                  **kwargs
                  ):
 
@@ -904,6 +932,7 @@ def __init__(self, ax,
             filternorm=filternorm,
             filterrad=filterrad,
             resample=resample,
+            interpolation_stage=interpolation_stage,
             **kwargs
         )
 

lib/matplotlib/pyplot.py

@@ -2864,12 +2864,13 @@ def hlines(
 def imshow(
         X, cmap=None, norm=None, aspect=None, interpolation=None,
         alpha=None, vmin=None, vmax=None, origin=None, extent=None, *,
-        filternorm=True, filterrad=4.0, resample=None, url=None,
-        data=None, **kwargs):
+        interpolation_stage=None, filternorm=True, filterrad=4.0,
+        resample=None, url=None, data=None, **kwargs):
     __ret = gca().imshow(
         X, cmap=cmap, norm=norm, aspect=aspect,
         interpolation=interpolation, alpha=alpha, vmin=vmin,
         vmax=vmax, origin=origin, extent=extent,
+        interpolation_stage=interpolation_stage,
         filternorm=filternorm, filterrad=filterrad, resample=resample,
         url=url, **({"data": data} if data is not None else {}),
         **kwargs)

lib/matplotlib/tests/baseline_images/test_image/rgba_antialias.png

@@ -1328,3 +1328,58 @@ def test_nonuniform_and_pcolor():
         ax.set_axis_off()
         # NonUniformImage "leaks" out of extents, not PColorImage.
         ax.set(xlim=(0, 10))
+
+
+@image_comparison(["rgba_antialias.png"], style="mpl20",
+                  remove_text=True)
+def test_rgba_antialias():
+    fig, axs = plt.subplots(2, 2, figsize=(3.5, 3.5), sharex=False,
+                            sharey=False, constrained_layout=True)
+    N = 250
+    aa = np.ones((N, N))
+    aa[::2, :] = -1
+
+    x = np.arange(N) / N - 0.5
+    y = np.arange(N) / N - 0.5
+
+    X, Y = np.meshgrid(x, y)
+    R = np.sqrt(X**2 + Y**2)
+    f0 = 10
+    k = 75
+    # aliased concentric circles
+    a = np.sin(np.pi * 2 * (f0 * R + k * R**2 / 2))
+
+    # stripes on lhs
+    a[:int(N/2), :][R[:int(N/2), :] < 0.4] = -1
+    a[:int(N/2), :][R[:int(N/2), :] < 0.3] = 1
+    aa[:, int(N/2):] = a[:, int(N/2):]
+
+    # set some over/unders and NaNs
+    aa[20:50, 20:50] = np.NaN
+    aa[70:90, 70:90] = 1e6
+    aa[70:90, 20:30] = -1e6
+    aa[70:90, 195:215] = 1e6
+    aa[20:30, 195:215] = -1e6
+
+    cmap = copy(plt.cm.RdBu_r)
+    cmap.set_over('yellow')
+    cmap.set_under('cyan')
+
+    axs = axs.flatten()
+    # zoom in
+    axs[0].imshow(aa, interpolation='nearest', cmap=cmap, vmin=-1.2, vmax=1.2)
+    axs[0].set_xlim([N/2-25, N/2+25])
+    axs[0].set_ylim([N/2+50, N/2-10])
+
+    # no anti-alias
+    axs[1].imshow(aa, interpolation='nearest', cmap=cmap, vmin=-1.2, vmax=1.2)
+
+    # data antialias: Note no purples, and white in circle.  Note
+    # that alternating red and blue stripes become white.
+    axs[2].imshow(aa, interpolation='antialiased', interpolation_stage='data',
+                  cmap=cmap, vmin=-1.2, vmax=1.2)
+
+    # rgba antialias: Note purples at boundary with circle.  Note that
+    # alternating red and blue stripes become purple
+    axs[3].imshow(aa, interpolation='antialiased', interpolation_stage='rgba',
+                  cmap=cmap, vmin=-1.2, vmax=1.2)