diff --git a/examples/images_contours_and_fields/image_transparency_blend.py b/examples/images_contours_and_fields/image_transparency_blend.py
@@ -6,12 +6,10 @@
 Blend transparency with color to highlight parts of data with imshow.
 
 A common use for :func:`matplotlib.pyplot.imshow` is to plot a 2-D statistical
-map. While ``imshow`` makes it easy to visualize a 2-D matrix as an image,
-it doesn't easily let you add transparency to the output. For example, one can
-plot a statistic (such as a t-statistic) and color the transparency of
-each pixel according to its p-value. This example demonstrates how you can
-achieve this effect using :class:`matplotlib.colors.Normalize`. Note that it is
-not possible to directly pass alpha values to :func:`matplotlib.pyplot.imshow`.
+map. The function makes it easy to visualize a 2-D matrix as an image and add
+transparency to the output. For example, one can plot a statistic (such as a
+t-statistic) and color the transparency of each pixel according to its p-value.
+This example demonstrates how you can achieve this effect.
 
 First we will generate some data, in this case, we'll create two 2-D "blobs"
 in a 2-D grid. One blob will be positive, and the other negative.
@@ -50,42 +48,38 @@ def normal_pdf(x, mean, var):
 # We'll also create a grey background into which the pixels will fade
 greys = np.full((*weights.shape, 3), 70, dtype=np.uint8)
 
-# First we'll plot these blobs using only ``imshow``.
+# First we'll plot these blobs using ``imshow`` without transparency.
 vmax = np.abs(weights).max()
-vmin = -vmax
-cmap = plt.cm.RdYlBu
+imshow_kwargs = {
+    'vmax': vmax,
+    'vmin': -vmax,
+    'cmap': 'RdYlBu',
+    'extent': (xmin, xmax, ymin, ymax),
+}
 
 fig, ax = plt.subplots()
 ax.imshow(greys)
-ax.imshow(weights, extent=(xmin, xmax, ymin, ymax), cmap=cmap)
+ax.imshow(weights, **imshow_kwargs)
 ax.set_axis_off()
 
 ###############################################################################
 # Blending in transparency
 # ========================
 #
 # The simplest way to include transparency when plotting data with
-# :func:`matplotlib.pyplot.imshow` is to convert the 2-D data array to a
-# 3-D image array of rgba values. This can be done with
-# :class:`matplotlib.colors.Normalize`. For example, we'll create a gradient
+# :func:`matplotlib.pyplot.imshow` is to pass an array matching the shape of
+# the data to the ``alpha`` argument. For example, we'll create a gradient
 # moving from left to right below.
 
 # Create an alpha channel of linearly increasing values moving to the right.
 alphas = np.ones(weights.shape)
 alphas[:, 30:] = np.linspace(1, 0, 70)
 
-# Normalize the colors b/w 0 and 1, we'll then pass an MxNx4 array to imshow
-colors = Normalize(vmin, vmax, clip=True)(weights)
-colors = cmap(colors)
-
-# Now set the alpha channel to the one we created above
-colors[..., -1] = alphas
-
 # Create the figure and image
 # Note that the absolute values may be slightly different
 fig, ax = plt.subplots()
 ax.imshow(greys)
-ax.imshow(colors, extent=(xmin, xmax, ymin, ymax))
+ax.imshow(weights, alpha=alphas, **imshow_kwargs)
 ax.set_axis_off()
 
 ###############################################################################
@@ -102,18 +96,11 @@ def normal_pdf(x, mean, var):
 alphas = Normalize(0, .3, clip=True)(np.abs(weights))
 alphas = np.clip(alphas, .4, 1)  # alpha value clipped at the bottom at .4
 
-# Normalize the colors b/w 0 and 1, we'll then pass an MxNx4 array to imshow
-colors = Normalize(vmin, vmax)(weights)
-colors = cmap(colors)
-
-# Now set the alpha channel to the one we created above
-colors[..., -1] = alphas
-
 # Create the figure and image
 # Note that the absolute values may be slightly different
 fig, ax = plt.subplots()
 ax.imshow(greys)
-ax.imshow(colors, extent=(xmin, xmax, ymin, ymax))
+ax.imshow(weights, alpha=alphas, **imshow_kwargs)
 
 # Add contour lines to further highlight different levels.
 ax.contour(weights[::-1], levels=[-.1, .1], colors='k', linestyles='-')

diff --git a/lib/matplotlib/axes/_axes.py b/lib/matplotlib/axes/_axes.py
@@ -5504,9 +5504,10 @@ 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*.
 
-        alpha : scalar, optional
+        alpha : scalar or array-like, optional
             The alpha blending value, between 0 (transparent) and 1 (opaque).
-            This parameter is ignored for RGBA input data.
+            If *alpha* is an array, the alpha blending values are applied pixel
+            by pixel, and *alpha* must have the same shape as *X*.
 
         vmin, vmax : scalar, optional
             When using scalar data and no explicit *norm*, *vmin* and *vmax*

diff --git a/lib/matplotlib/image.py b/lib/matplotlib/image.py
@@ -7,6 +7,7 @@
 import math
 import os
 import logging
+from numbers import Number
 from pathlib import Path
 import urllib.parse
 
@@ -95,7 +96,7 @@ def composite_images(images, renderer, magnification=1.0):
         if data is not None:
             x *= magnification
             y *= magnification
-            parts.append((data, x, y, image.get_alpha() or 1.0))
+            parts.append((data, x, y, image._get_scalar_alpha()))
             bboxes.append(
                 Bbox([[x, y], [x + data.shape[1], y + data.shape[0]]]))
 
@@ -281,9 +282,29 @@ def set_alpha(self, alpha):
         ----------
         alpha : float
         """
-        martist.Artist.set_alpha(self, alpha)
+        if alpha is not None and not isinstance(alpha, Number):
+            alpha = np.asarray(alpha)
+            if alpha.ndim != 2:
+                raise TypeError('alpha must be a float, two-dimensional '
+                                'array, or None')
+        self._alpha = alpha
+        self.pchanged()
+        self.stale = True
         self._imcache = None
 
+    def _get_scalar_alpha(self):
+        """
+        Get a scalar alpha value to be applied to the artist as a whole.
+
+        If the alpha value is a matrix, the method returns 1.0 because pixels
+        have individual alpha values (see `~._ImageBase._make_image` for
+        details). If the alpha value is a scalar, the method returns said value
+        to be applied to the artist as a whole because pixels do not have
+        individual alpha values.
+        """
+        return 1.0 if self._alpha is None or np.ndim(self._alpha) > 0 \
+            else self._alpha
+
     def changed(self):
         """
         Call this whenever the mappable is changed so observers can
@@ -487,14 +508,17 @@ def _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification=1.0,
                 # pixel it will be between [0, 1] (such as a rotated image).
                 out_mask = np.isnan(out_alpha)
                 out_alpha[out_mask] = 1
+                # Apply the pixel-by-pixel alpha values if present
+                alpha = self.get_alpha()
+                if alpha is not None and np.ndim(alpha) > 0:
+                    out_alpha *= _resample(self, alpha, out_shape,
+                                           t, resample=True)
                 # mask and run through the norm
                 output = self.norm(np.ma.masked_array(A_resampled, out_mask))
             else:
                 if A.shape[2] == 3:
                     A = _rgb_to_rgba(A)
-                alpha = self.get_alpha()
-                if alpha is None:
-                    alpha = 1
+                alpha = self._get_scalar_alpha()
                 output_alpha = _resample(  # resample alpha channel
                     self, A[..., 3], out_shape, t, alpha=alpha)
                 output = _resample(  # resample rgb channels
@@ -509,9 +533,7 @@ def _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification=1.0,
 
             # Apply alpha *after* if the input was greyscale without a mask
             if A.ndim == 2:
-                alpha = self.get_alpha()
-                if alpha is None:
-                    alpha = 1
+                alpha = self._get_scalar_alpha()
                 alpha_channel = output[:, :, 3]
                 alpha_channel[:] = np.asarray(
                     np.asarray(alpha_channel, np.float32) * out_alpha * alpha,
@@ -593,7 +615,7 @@ def draw(self, renderer, *args, **kwargs):
         # actually render the image.
         gc = renderer.new_gc()
         self._set_gc_clip(gc)
-        gc.set_alpha(self.get_alpha())
+        gc.set_alpha(self._get_scalar_alpha())
         gc.set_url(self.get_url())
         gc.set_gid(self.get_gid())
 

diff --git a/lib/matplotlib/tests/test_image.py b/lib/matplotlib/tests/test_image.py
@@ -1118,3 +1118,22 @@ def test_image_cursor_formatting():
 
     data = np.nan
     assert im.format_cursor_data(data) == '[nan]'
+
+
+@check_figures_equal()
+def test_image_array_alpha(fig_test, fig_ref):
+    '''per-pixel alpha channel test'''
+    x = np.linspace(0, 1)
+    xx, yy = np.meshgrid(x, x)
+
+    zz = np.exp(- 3 * ((xx - 0.5) ** 2) + (yy - 0.7 ** 2))
+    alpha = zz / zz.max()
+
+    cmap = plt.get_cmap('viridis')
+    ax = fig_test.add_subplot(111)
+    ax.imshow(zz, alpha=alpha, cmap=cmap, interpolation='nearest')
+
+    ax = fig_ref.add_subplot(111)
+    rgba = cmap(colors.Normalize()(zz))
+    rgba[..., -1] = alpha
+    ax.imshow(rgba, interpolation='nearest')