ENH: anti-alias down-sampled images

Author: jklymak

.flake8

@@ -94,6 +94,7 @@ per-file-ignores =
     examples/images_contours_and_fields/contourf_hatching.py: E402
     examples/images_contours_and_fields/contourf_log.py: E402
     examples/images_contours_and_fields/demo_bboximage.py: E402
+    examples/images_contours_and_fields/image_antialiasing.py: E402
     examples/images_contours_and_fields/image_clip_path.py: E402
     examples/images_contours_and_fields/image_demo.py: E402
     examples/images_contours_and_fields/image_masked.py: E402

doc/api/next_api_changes/2019-07-17-JMK.rst

@@ -0,0 +1,22 @@
+Default interpolation for `image` is new "antialiased" option
+-------------------------------------------------------------
+
+Images displayed in Matplotlib previously used nearest-neighbor
+interpolation, leading to aliasing effects for downscaling and non-integer
+upscaling.  
+
+New default for :rc:`image.interpolation` is the new option "antialiased".
+`imshow(A, interpolation='antialiased')` will apply a Hanning filter when
+resampling the data in A for display (or saving to file) *if* the upsample
+rate is less than a factor of three, and not an integer; downsampled data is
+always smoothed at resampling.
+
+To get the old behavior, set :rc:`interpolation` to the old default "nearest"
+(or specify the ``interpolation`` kwarg of `.Axes.imshow`)
+
+To always get the anti-aliasing behavior, no matter what the up/down sample
+rate, set :rc:`interpolation` to "hanning" (or one of the other filters
+available.
+
+Note that the "hanning" filter was chosen because it has only a modest
+performance penalty.  Anti-aliasing can be improved with other filters.

examples/images_contours_and_fields/image_antialiasing.py

@@ -0,0 +1,81 @@
+"""
+==================
+Image Antialiasing
+==================
+
+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.
+
+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.
+
+Other anti-aliasing filters can be specified in `.Axes.imshow` using the
+*interpolation* kwarg.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+###############################################################################
+# First we generate an image with varying frequency content:
+x = np.arange(500) / 500 - 0.5
+y = np.arange(500) / 500 - 0.5
+
+X, Y = np.meshgrid(x, y)
+R = np.sqrt(X**2 + Y**2)
+f0 = 10
+k = 250
+a = np.sin(np.pi * 2 * (f0 * R + k * R**2 / 2))
+
+
+###############################################################################
+# The following images are subsampled from 1000 data pixels to 604 rendered
+# pixels. The Moire patterns in the "nearest" interpolation are caused by the
+# high-frequency data being subsampled.  The "antialiased" image
+# 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 n, interp in enumerate(['nearest', 'antialiased']):
+    im = axs[n].imshow(a, interpolation=interp, cmap='gray')
+    axs[n].set_title(interp)
+plt.show()
+
+###############################################################################
+# Even up-sampling an image will lead to Moire patterns unless the upsample
+# is an integer number of pixels.
+fig, ax = plt.subplots(1, 1, figsize=(5.3, 5.3))
+ax.set_position([0, 0, 1, 1])
+im = ax.imshow(a, interpolation='nearest', cmap='gray')
+plt.show()
+
+###############################################################################
+# The patterns aren't as bad, but still benefit from anti-aliasing
+fig, ax = plt.subplots(1, 1, figsize=(5.3, 5.3))
+ax.set_position([0, 0, 1, 1])
+im = ax.imshow(a, interpolation='antialiased', cmap='gray')
+plt.show()
+
+###############################################################################
+# If the small Moire patterns in the default "hanning" antialiasing are
+# still undesireable, then we can use other filters.
+fig, axs = plt.subplots(1, 2, figsize=(7, 4), constrained_layout=True)
+for n, interp in enumerate(['hanning', 'lanczos']):
+    im = axs[n].imshow(a, interpolation=interp, cmap='gray')
+    axs[n].set_title(interp)
+plt.show()
+
+
+#############################################################################
+#
+# ------------
+#
+# References
+# """"""""""
+#
+# The use of the following functions and methods is shown
+# in this example:
+
+import matplotlib
+matplotlib.axes.Axes.imshow

examples/images_contours_and_fields/interpolation_methods.py

@@ -9,11 +9,14 @@
 If `interpolation` is None, it defaults to the :rc:`image.interpolation`
 (default: ``'nearest'``). If the interpolation is ``'none'``, then no
 interpolation is performed for the Agg, ps and pdf backends. Other backends
-will default to ``'nearest'``.
+will default to ``'antialiased'``.
 
 For the Agg, ps and pdf backends, ``interpolation = 'none'`` works well when a
 big image is scaled down, while ``interpolation = 'nearest'`` works well when
 a small image is scaled up.
+
+See :doc:`/gallery/images_contours_and_fields/image_antialiasing` for a
+discussion on the default `interpolation="antialiased"` option.
 """
 
 import matplotlib.pyplot as plt

lib/matplotlib/axes/_axes.py

@@ -5488,20 +5488,30 @@ def imshow(self, X, cmap=None, norm=None, aspect=None,
             The interpolation method used. If *None*
             :rc:`image.interpolation` is used, which defaults to 'nearest'.
 
-            Supported values are 'none', 'nearest', 'bilinear', 'bicubic',
-            'spline16', 'spline36', 'hanning', 'hamming', 'hermite', 'kaiser',
-            'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell', 'sinc',
-            'lanczos'.
+            Supported values are 'none', 'antialiased', 'nearest', 'bilinear',
+            'bicubic', 'spline16', 'spline36', 'hanning', 'hamming', 'hermite',
+            'kaiser', 'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell',
+            'sinc', 'lanczos'.
 
             If *interpolation* is 'none', then no interpolation is performed
             on the Agg, ps, pdf and svg backends. Other backends will fall back
             to 'nearest'. Note that most SVG renders perform interpolation at
             rendering and that the default interpolation method they implement
             may differ.
 
+            If *interpolation* is the default 'antialiased', then 'nearest'
+            interpolation is used if the image is upsampled by more than a
+            factor of three (i.e. the number of display pixels is at least
+            three times the size of the data array).  If the upsampling rate is
+            smaller than 3, or the image is downsampled, then 'hanning'
+            interpolation is used to act as an anti-aliasing filter, unless the
+            image happens to be upsampled by exactly a factor of two or one.
+
             See
             :doc:`/gallery/images_contours_and_fields/interpolation_methods`
-            for an overview of the supported interpolation methods.
+            for an overview of the supported interpolation methods, and
+            :doc:`/gallery/images_contours_and_fields/image_antialiasing` for
+            a discussion of image antialiasing.
 
             Some interpolation methods require an additional radius parameter,
             which can be set by *filterrad*. Additionally, the antigrain image

lib/matplotlib/image.py

@@ -33,6 +33,7 @@
 
 # map interpolation strings to module constants
 _interpd_ = {
+    'antialiased': _image.NEAREST,  # this will use nearest or Hanning...
     'none': _image.NEAREST,  # fall back to nearest when not supported
     'nearest': _image.NEAREST,
     'bilinear': _image.BILINEAR,
@@ -168,11 +169,34 @@ def _resample(
     allocating the output array and fetching the relevant properties from the
     Image object *image_obj*.
     """
+
+    # decide if we need to apply anti-aliasing if the data is upsampled:
+    # compare the number of displayed pixels to the number of
+    # the data pixels.
+    interpolation = image_obj.get_interpolation()
+    if interpolation == 'antialiased':
+        # don't antialias if upsampling by an integer number or
+        # if zooming in more than a factor of 3
+        shape = list(data.shape)
+        if image_obj.origin == 'upper':
+            shape[0] = 0
+        dispx, dispy = transform.transform([shape[1], shape[0]])
+
+        if ((dispx > 3 * data.shape[1] or
+                dispx == data.shape[1] or
+                dispx == 2 * data.shape[1]) and
+            (dispy > 3 * data.shape[0] or
+                dispy == data.shape[0] or
+                dispy == 2 * data.shape[0])):
+            interpolation = 'nearest'
+        else:
+            interpolation = 'hanning'
+
     out = np.zeros(out_shape + data.shape[2:], data.dtype)  # 2D->2D, 3D->3D.
     if resample is None:
         resample = image_obj.get_resample()
     _image.resample(data, out, transform,
-                    _interpd_[image_obj.get_interpolation()],
+                    _interpd_[interpolation],
                     resample,
                     alpha,
                     image_obj.get_filternorm(),
@@ -432,7 +456,6 @@ def _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification=1.0,
                 A_scaled += 0.1
                 # resample the input data to the correct resolution and shape
                 A_resampled = _resample(self, A_scaled, out_shape, t)
-
                 # done with A_scaled now, remove from namespace to be sure!
                 del A_scaled
                 # un-scale the resampled data to approximately the
@@ -690,9 +713,10 @@ def get_interpolation(self):
         """
         Return the interpolation method the image uses when resizing.
 
-        One of 'nearest', 'bilinear', 'bicubic', 'spline16', 'spline36',
-        'hanning', 'hamming', 'hermite', 'kaiser', 'quadric', 'catrom',
-        'gaussian', 'bessel', 'mitchell', 'sinc', 'lanczos', or 'none'.
+        One of 'antialiased', 'nearest', 'bilinear', 'bicubic', 'spline16',
+        'spline36', 'hanning', 'hamming', 'hermite', 'kaiser', 'quadric',
+        'catrom', 'gaussian', 'bessel', 'mitchell', 'sinc', 'lanczos',
+        or 'none'.
 
         """
         return self._interpolation
@@ -708,9 +732,9 @@ def set_interpolation(self, s):
 
         Parameters
         ----------
-        s : {'nearest', 'bilinear', 'bicubic', 'spline16', 'spline36', \
-'hanning', 'hamming', 'hermite', 'kaiser', 'quadric', 'catrom', 'gaussian', \
-'bessel', 'mitchell', 'sinc', 'lanczos', 'none'}
+        s : {'antialiased', 'nearest', 'bilinear', 'bicubic', 'spline16',
+'spline36', 'hanning', 'hamming', 'hermite', 'kaiser', 'quadric', 'catrom', \
+'gaussian', 'bessel', 'mitchell', 'sinc', 'lanczos', 'none'}
 
         """
         if s is None:

lib/matplotlib/rcsetup.py

@@ -1100,7 +1100,7 @@ def _validate_linestyle(ls):
     'mathtext.fallback_to_cm': [True, validate_bool],
 
     'image.aspect':        ['equal', validate_aspect],  # equal, auto, a number
-    'image.interpolation': ['nearest', validate_string],
+    'image.interpolation': ['antialiased', validate_string],
     'image.cmap':          ['viridis', validate_string],  # gray, jet, etc.
     'image.lut':           [256, validate_int],  # lookup table
     'image.origin':        ['upper',

lib/matplotlib/tests/test_axes.py

@@ -915,6 +915,8 @@ def test_nonfinite_limits():
 
 @image_comparison(['imshow', 'imshow'], remove_text=True, style='mpl20')
 def test_imshow():
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
     # Create a NxN image
     N = 100
     (x, y) = np.indices((N, N))
@@ -936,6 +938,8 @@ def test_imshow():
 @image_comparison(['imshow_clip'], style='mpl20')
 def test_imshow_clip():
     # As originally reported by Gellule Xg <gellule.xg@free.fr>
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
 
     # Create a NxN image
     N = 100
@@ -3732,6 +3736,10 @@ def test_subplot_key_hash():
                   remove_text=True, tol=0.07, style='default')
 def test_specgram_freqs():
     '''test axes.specgram in default (psd) mode with sinusoidal stimuli'''
+
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
+
     n = 1000
     Fs = 10.
 
@@ -3784,6 +3792,10 @@ def test_specgram_freqs():
                   remove_text=True, tol=0.01, style='default')
 def test_specgram_noise():
     '''test axes.specgram in default (psd) mode with noise stimuli'''
+
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
+
     np.random.seed(0)
 
     n = 1000
@@ -3831,6 +3843,10 @@ def test_specgram_noise():
                   remove_text=True, tol=0.07, style='default')
 def test_specgram_magnitude_freqs():
     '''test axes.specgram in magnitude mode with sinusoidal stimuli'''
+
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
+
     n = 1000
     Fs = 10.
 
@@ -3886,6 +3902,10 @@ def test_specgram_magnitude_freqs():
                   remove_text=True, style='default')
 def test_specgram_magnitude_noise():
     '''test axes.specgram in magnitude mode with noise stimuli'''
+
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
+
     np.random.seed(0)
 
     n = 1000
@@ -3932,6 +3952,10 @@ def test_specgram_magnitude_noise():
                   remove_text=True, tol=0.007, style='default')
 def test_specgram_angle_freqs():
     '''test axes.specgram in angle mode with sinusoidal stimuli'''
+
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
+
     n = 1000
     Fs = 10.
 
@@ -3986,6 +4010,10 @@ def test_specgram_angle_freqs():
                   remove_text=True, style='default')
 def test_specgram_noise_angle():
     '''test axes.specgram in angle mode with noise stimuli'''
+
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
+
     np.random.seed(0)
 
     n = 1000
@@ -4032,6 +4060,9 @@ def test_specgram_noise_angle():
                   remove_text=True, style='default')
 def test_specgram_freqs_phase():
     '''test axes.specgram in phase mode with sinusoidal stimuli'''
+
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
     n = 1000
     Fs = 10.
 
@@ -4086,6 +4117,9 @@ def test_specgram_freqs_phase():
                   remove_text=True, style='default')
 def test_specgram_noise_phase():
     '''test axes.specgram in phase mode with noise stimuli'''
+
+    # use former defaults to match existing baseline image
+    matplotlib.rcParams['image.interpolation'] = 'nearest'
     np.random.seed(0)
 
     n = 1000