diff --git a/doc/users/whats_new/axes3d_voxels.rst b/doc/users/whats_new/axes3d_voxels.rst
new file mode 100644
index 000000000000..2f22d77b81bf
--- /dev/null
+++ b/doc/users/whats_new/axes3d_voxels.rst
@@ -0,0 +1,5 @@
+``voxels`` function for mplot3d
+-------------------------------
+:class:`~mpl_toolkits.mplot3d.axes3d.Axes3D` now has a ``voxels`` method, for
+visualizing boolean 3d data. Uses could include plotting a sparse 3D heat map,
+or visualizing a volumetric model.
diff --git a/examples/mplot3d/voxels.py b/examples/mplot3d/voxels.py
new file mode 100644
index 000000000000..76cf64c33a00
--- /dev/null
+++ b/examples/mplot3d/voxels.py
@@ -0,0 +1,35 @@
+'''
+==========================
+3D voxel / volumetric plot
+==========================
+
+Demonstrates plotting 3D volumetric objects with ``ax.voxels``
+'''
+
+import matplotlib.pyplot as plt
+import numpy as np
+from mpl_toolkits.mplot3d import Axes3D
+
+# prepare some coordinates
+x, y, z = np.indices((8, 8, 8))
+
+# draw cuboids in the top left and bottom right corners, and a link between them
+cube1 = (x < 3) & (y < 3) & (z < 3)
+cube2 = (x >= 5) & (y >= 5) & (z >= 5)
+link = abs(x - y) + abs(y - z) + abs(z - x) <= 2
+
+# combine the objects into a single boolean array
+voxels = cube1 | cube2 | link
+
+# set the colors of each object
+colors = np.empty(voxels.shape, dtype=object)
+colors[link] = 'red'
+colors[cube1] = 'blue'
+colors[cube2] = 'green'
+
+# and plot everything
+fig = plt.figure()
+ax = fig.gca(projection='3d')
+ax.voxels(voxels, facecolors=colors, edgecolor='k')
+
+plt.show()
diff --git a/examples/mplot3d/voxels_numpy_logo.py b/examples/mplot3d/voxels_numpy_logo.py
new file mode 100644
index 000000000000..648a3cff7822
--- /dev/null
+++ b/examples/mplot3d/voxels_numpy_logo.py
@@ -0,0 +1,47 @@
+'''
+===============================
+3D voxel plot of the numpy logo
+===============================
+
+Demonstrates using ``ax.voxels`` with uneven coordinates
+'''
+import matplotlib.pyplot as plt
+import numpy as np
+from mpl_toolkits.mplot3d import Axes3D
+
+
+def explode(data):
+    size = np.array(data.shape)*2
+    data_e = np.zeros(size - 1, dtype=data.dtype)
+    data_e[::2, ::2, ::2] = data
+    return data_e
+
+# build up the numpy logo
+n_voxels = np.zeros((4, 3, 4), dtype=bool)
+n_voxels[0, 0, :] = True
+n_voxels[-1, 0, :] = True
+n_voxels[1, 0, 2] = True
+n_voxels[2, 0, 1] = True
+facecolors = np.where(n_voxels, '#FFD65DC0', '#7A88CCC0')
+edgecolors = np.where(n_voxels, '#BFAB6E', '#7D84A6')
+filled = np.ones(n_voxels.shape)
+
+# upscale the above voxel image, leaving gaps
+filled_2 = explode(filled)
+fcolors_2 = explode(facecolors)
+ecolors_2 = explode(edgecolors)
+
+# Shrink the gaps
+x, y, z = np.indices(np.array(filled_2.shape) + 1).astype(float) // 2
+x[0::2, :, :] += 0.05
+y[:, 0::2, :] += 0.05
+z[:, :, 0::2] += 0.05
+x[1::2, :, :] += 0.95
+y[:, 1::2, :] += 0.95
+z[:, :, 1::2] += 0.95
+
+fig = plt.figure()
+ax = fig.gca(projection='3d')
+ax.voxels(x, y, z, filled_2, facecolors=fcolors_2, edgecolors=ecolors_2)
+
+plt.show()
diff --git a/examples/mplot3d/voxels_rgb.py b/examples/mplot3d/voxels_rgb.py
new file mode 100644
index 000000000000..1b577cad47fe
--- /dev/null
+++ b/examples/mplot3d/voxels_rgb.py
@@ -0,0 +1,45 @@
+'''
+==========================================
+3D voxel / volumetric plot with rgb colors
+==========================================
+
+Demonstrates using ``ax.voxels`` to visualize parts of a color space
+'''
+
+import matplotlib.pyplot as plt
+import numpy as np
+from mpl_toolkits.mplot3d import Axes3D
+
+
+def midpoints(x):
+    sl = ()
+    for i in range(x.ndim):
+        x = (x[sl + np.index_exp[:-1]] + x[sl + np.index_exp[1:]]) / 2.0
+        sl += np.index_exp[:]
+    return x
+
+# prepare some coordinates, and attach rgb values to each
+r, g, b = np.indices((17, 17, 17)) / 16.0
+rc = midpoints(r)
+gc = midpoints(g)
+bc = midpoints(b)
+
+# define a sphere about [0.5, 0.5, 0.5]
+sphere = (rc - 0.5)**2 + (gc - 0.5)**2 + (bc - 0.5)**2 < 0.5**2
+
+# combine the color components
+colors = np.zeros(sphere.shape + (3,))
+colors[..., 0] = rc
+colors[..., 1] = gc
+colors[..., 2] = bc
+
+# and plot everything
+fig = plt.figure()
+ax = fig.gca(projection='3d')
+ax.voxels(r, g, b, sphere,
+          facecolors=colors,
+          edgecolors=np.clip(2*colors - 0.5, 0, 1),  # brighter
+          linewidth=0.5)
+ax.set(xlabel='r', ylabel='g', zlabel='b')
+
+plt.show()
diff --git a/examples/mplot3d/voxels_torus.py b/examples/mplot3d/voxels_torus.py
new file mode 100644
index 000000000000..4f60e31403d8
--- /dev/null
+++ b/examples/mplot3d/voxels_torus.py
@@ -0,0 +1,47 @@
+'''
+=======================================================
+3D voxel / volumetric plot with cylindrical coordinates
+=======================================================
+
+Demonstrates using the ``x, y, z`` arguments of ``ax.voxels``.
+'''
+
+import matplotlib.pyplot as plt
+import matplotlib.colors
+import numpy as np
+from mpl_toolkits.mplot3d import Axes3D
+
+
+def midpoints(x):
+    sl = ()
+    for i in range(x.ndim):
+        x = (x[sl + np.index_exp[:-1]] + x[sl + np.index_exp[1:]]) / 2.0
+        sl += np.index_exp[:]
+    return x
+
+# prepare some coordinates, and attach rgb values to each
+r, theta, z = np.mgrid[0:1:11j, 0:np.pi*2:25j, -0.5:0.5:11j]
+x = r*np.cos(theta)
+y = r*np.sin(theta)
+
+rc, thetac, zc = midpoints(r), midpoints(theta), midpoints(z)
+
+# define a wobbly torus about [0.7, *, 0]
+sphere = (rc - 0.7)**2 + (zc + 0.2*np.cos(thetac*2))**2 < 0.2**2
+
+# combine the color components
+hsv = np.zeros(sphere.shape + (3,))
+hsv[..., 0] = thetac / (np.pi*2)
+hsv[..., 1] = rc
+hsv[..., 2] = zc + 0.5
+colors = matplotlib.colors.hsv_to_rgb(hsv)
+
+# and plot everything
+fig = plt.figure()
+ax = fig.gca(projection='3d')
+ax.voxels(x, y, z, sphere,
+          facecolors=colors,
+          edgecolors=np.clip(2*colors - 0.5, 0, 1),  # brighter
+          linewidth=0.5)
+
+plt.show()
diff --git a/lib/mpl_toolkits/mplot3d/axes3d.py b/lib/mpl_toolkits/mplot3d/axes3d.py
index e1e968430f5f..52d8ed842966 100644
--- a/lib/mpl_toolkits/mplot3d/axes3d.py
+++ b/lib/mpl_toolkits/mplot3d/axes3d.py
@@ -17,6 +17,7 @@
 from six.moves import map, xrange, zip, reduce
 
 import warnings
+from collections import defaultdict
 
 import numpy as np
 import matplotlib.axes as maxes
@@ -30,6 +31,7 @@
 from matplotlib.tri.triangulation import Triangulation
 from matplotlib import colors as mcolors
 from matplotlib.colors import Normalize, LightSource
+from matplotlib.cbook._backports import broadcast_to
 
 from . import art3d
 from . import proj3d
@@ -2744,6 +2746,212 @@ def calc_arrow(uvw, angle=15):
 
     quiver3D = quiver
 
+    def voxels(self, *args, **kwargs):
+        """
+        ax.voxels([x, y, z,] /, filled, **kwargs)
+
+        Plot a set of filled voxels
+
+        All voxels are plotted as 1x1x1 cubes on the axis, with filled[0,0,0]
+        placed with its lower corner at the origin. Occluded faces are not
+        plotted.
+
+        Call signatures::
+
+            voxels(filled, facecolors=fc, edgecolors=ec, **kwargs)
+            voxels(x, y, z, filled, facecolors=fc, edgecolors=ec, **kwargs)
+
+        .. versionadded:: 2.1
+
+        Parameters
+        ----------
+        filled : 3D np.array of bool
+            A 3d array of values, with truthy values indicating which voxels
+            to fill
+
+        x, y, z : 3D np.array, optional
+            The coordinates of the corners of the voxels. This should broadcast
+            to a shape one larger in every dimension than the shape of `filled`.
+            These can be used to plot non-cubic voxels.
+
+            If not specified, defaults to increasing integers along each axis,
+            like those returned by :func:`~numpy.indices`.
+            As indicated by the ``/`` in the function signature, these arguments
+            can only be passed positionally.
+
+        facecolors, edgecolors : array_like, optional
+            The color to draw the faces and edges of the voxels. Can only be
+            passed as keyword arguments.
+            This parameter can be:
+
+              - A single color value, to color all voxels the same color. This
+                can be either a string, or a 1D rgb/rgba array
+              - ``None``, the default, to use a single color for the faces, and
+                the style default for the edges.
+              - A 3D ndarray of color names, with each item the color for the
+                corresponding voxel. The size must match the voxels.
+              - A 4D ndarray of rgb/rgba data, with the components along the
+                last axis.
+
+        **kwargs
+            Additional keyword arguments to pass onto
+            :func:`~mpl_toolkits.mplot3d.art3d.Poly3DCollection`
+
+        Returns
+        -------
+        faces : dict
+            A dictionary indexed by coordinate, where ``faces[i,j,k]`` is a
+            `Poly3DCollection` of the faces drawn for the voxel
+            ``filled[i,j,k]``. If no faces were drawn for a given voxel, either
+            because it was not asked to be drawn, or it is fully occluded, then
+            ``(i,j,k) not in faces``.
+
+        Examples
+        --------
+        .. plot:: gallery/mplot3d/voxels.py
+        .. plot:: gallery/mplot3d/voxels_rgb.py
+        .. plot:: gallery/mplot3d/voxels_torus.py
+        .. plot:: gallery/mplot3d/voxels_numpy_logo.py
+        """
+
+        # work out which signature we should be using, and use it to parse
+        # the arguments. Name must be voxels for the correct error message
+        if len(args) >= 3:
+            # underscores indicate position only
+            def voxels(__x, __y, __z, filled, **kwargs):
+                return (__x, __y, __z), filled, kwargs
+        else:
+            def voxels(filled, **kwargs):
+                return None, filled, kwargs
+
+        xyz, filled, kwargs = voxels(*args, **kwargs)
+
+        # check dimensions
+        if filled.ndim != 3:
+            raise ValueError("Argument filled must be 3-dimensional")
+        size = np.array(filled.shape, dtype=np.intp)
+
+        # check xyz coordinates, which are one larger than the filled shape
+        coord_shape = tuple(size + 1)
+        if xyz is None:
+            x, y, z = np.indices(coord_shape)
+        else:
+            x, y, z = (broadcast_to(c, coord_shape) for c in xyz)
+
+        def _broadcast_color_arg(color, name):
+            if np.ndim(color) in (0, 1):
+                # single color, like "red" or [1, 0, 0]
+                return broadcast_to(color, filled.shape + np.shape(color))
+            elif np.ndim(color) in (3, 4):
+                # 3D array of strings, or 4D array with last axis rgb
+                if np.shape(color)[:3] != filled.shape:
+                    raise ValueError(
+                        "When multidimensional, {} must match the shape of "
+                        "filled".format(name))
+                return color
+            else:
+                raise ValueError("Invalid {} argument".format(name))
+
+        # intercept the facecolors, handling defaults and broacasting
+        facecolors = kwargs.pop('facecolors', None)
+        if facecolors is None:
+            facecolors = self._get_patches_for_fill.get_next_color()
+        facecolors = _broadcast_color_arg(facecolors, 'facecolors')
+
+        # broadcast but no default on edgecolors
+        edgecolors = kwargs.pop('edgecolors', None)
+        edgecolors = _broadcast_color_arg(edgecolors, 'edgecolors')
+
+        # always scale to the full array, even if the data is only in the center
+        self.auto_scale_xyz(x, y, z)
+
+        # points lying on corners of a square
+        square = np.array([
+            [0, 0, 0],
+            [0, 1, 0],
+            [1, 1, 0],
+            [1, 0, 0]
+        ], dtype=np.intp)
+
+        voxel_faces = defaultdict(list)
+
+        def permutation_matrices(n):
+            """ Generator of cyclic permutation matices """
+            mat = np.eye(n, dtype=np.intp)
+            for i in range(n):
+                yield mat
+                mat = np.roll(mat, 1, axis=0)
+
+        # iterate over each of the YZ, ZX, and XY orientations, finding faces to
+        # render
+        for permute in permutation_matrices(3):
+            # find the set of ranges to iterate over
+            pc, qc, rc = permute.T.dot(size)
+            pinds = np.arange(pc)
+            qinds = np.arange(qc)
+            rinds = np.arange(rc)
+
+            square_rot = square.dot(permute.T)
+
+            # iterate within the current plane
+            for p in pinds:
+                for q in qinds:
+                    # iterate perpendicularly to the current plane, handling
+                    # boundaries. We only draw faces between a voxel and an
+                    # empty space, to avoid drawing internal faces.
+
+                    # draw lower faces
+                    p0 = permute.dot([p, q, 0])
+                    i0 = tuple(p0)
+                    if filled[i0]:
+                        voxel_faces[i0].append(p0 + square_rot)
+
+                    # draw middle faces
+                    for r1, r2 in zip(rinds[:-1], rinds[1:]):
+                        p1 = permute.dot([p, q, r1])
+                        p2 = permute.dot([p, q, r2])
+
+                        i1 = tuple(p1)
+                        i2 = tuple(p2)
+
+                        if filled[i1] and not filled[i2]:
+                            voxel_faces[i1].append(p2 + square_rot)
+                        elif not filled[i1] and filled[i2]:
+                            voxel_faces[i2].append(p2 + square_rot)
+
+                    # draw upper faces
+                    pk = permute.dot([p, q, rc-1])
+                    pk2 = permute.dot([p, q, rc])
+                    ik = tuple(pk)
+                    if filled[ik]:
+                        voxel_faces[ik].append(pk2 + square_rot)
+
+        # iterate over the faces, and generate a Poly3DCollection for each voxel
+        polygons = {}
+        for coord, faces_inds in voxel_faces.items():
+            # convert indices into 3D positions
+            if xyz is None:
+                faces = faces_inds
+            else:
+                faces = []
+                for face_inds in faces_inds:
+                    ind = face_inds[:, 0], face_inds[:, 1], face_inds[:, 2]
+                    face = np.empty(face_inds.shape)
+                    face[:, 0] = x[ind]
+                    face[:, 1] = y[ind]
+                    face[:, 2] = z[ind]
+                    faces.append(face)
+
+            poly = art3d.Poly3DCollection(faces,
+                facecolors=facecolors[coord],
+                edgecolors=edgecolors[coord],
+                **kwargs
+            )
+            self.add_collection3d(poly)
+            polygons[coord] = poly
+
+        return polygons
+
 
 def get_test_data(delta=0.05):
     '''
diff --git a/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-alpha.png b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-alpha.png
new file mode 100644
index 000000000000..ed84237b5450
Binary files /dev/null and b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-alpha.png differ
diff --git a/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-edge-style.png b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-edge-style.png
new file mode 100644
index 000000000000..c4dc28a988de
Binary files /dev/null and b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-edge-style.png differ
diff --git a/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-named-colors.png b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-named-colors.png
new file mode 100644
index 000000000000..6a139acc5e2d
Binary files /dev/null and b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-named-colors.png differ
diff --git a/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-rgb-data.png b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-rgb-data.png
new file mode 100644
index 000000000000..ab8308cead26
Binary files /dev/null and b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-rgb-data.png differ
diff --git a/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-simple.png b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-simple.png
new file mode 100644
index 000000000000..9889a775bfed
Binary files /dev/null and b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-simple.png differ
diff --git a/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-xyz.png b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-xyz.png
new file mode 100644
index 000000000000..4e2ecc61462e
Binary files /dev/null and b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/voxels-xyz.png differ
diff --git a/lib/mpl_toolkits/tests/test_mplot3d.py b/lib/mpl_toolkits/tests/test_mplot3d.py
index 9f74bd84223f..d5e65e81ce4d 100644
--- a/lib/mpl_toolkits/tests/test_mplot3d.py
+++ b/lib/mpl_toolkits/tests/test_mplot3d.py
@@ -569,6 +569,156 @@ def test_invalid_axes_limits(setter, side, value):
         getattr(obj, setter)(**limit)
 
 
+class TestVoxels(object):
+    @image_comparison(
+        baseline_images=['voxels-simple'],
+        extensions=['png'],
+        remove_text=True
+    )
+    def test_simple(self):
+        fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+        x, y, z = np.indices((5, 4, 3))
+        voxels = (x == y) | (y == z)
+        ax.voxels(voxels)
+
+    @image_comparison(
+        baseline_images=['voxels-edge-style'],
+        extensions=['png'],
+        remove_text=True,
+        style='default'
+    )
+    def test_edge_style(self):
+        fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+        x, y, z = np.indices((5, 5, 4))
+        voxels = ((x - 2)**2 + (y - 2)**2 + (z-1.5)**2) < 2.2**2
+        v = ax.voxels(voxels, linewidths=3, edgecolor='C1')
+
+        # change the edge color of one voxel
+        v[max(v.keys())].set_edgecolor('C2')
+
+    @image_comparison(
+        baseline_images=['voxels-named-colors'],
+        extensions=['png'],
+        remove_text=True
+    )
+    def test_named_colors(self):
+        """ test with colors set to a 3d object array of strings """
+        fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+        x, y, z = np.indices((10, 10, 10))
+        voxels = (x == y) | (y == z)
+        voxels = voxels & ~(x * y * z < 1)
+        colors = np.zeros((10, 10, 10), dtype=np.object_)
+        colors.fill('C0')
+        colors[(x < 5) & (y < 5)] = '0.25'
+        colors[(x + z) < 10] = 'cyan'
+        ax.voxels(voxels, facecolors=colors)
+
+    @image_comparison(
+        baseline_images=['voxels-rgb-data'],
+        extensions=['png'],
+        remove_text=True
+    )
+    def test_rgb_data(self):
+        """ test with colors set to a 4d float array of rgb data """
+        fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+        x, y, z = np.indices((10, 10, 10))
+        voxels = (x == y) | (y == z)
+        colors = np.zeros((10, 10, 10, 3))
+        colors[...,0] = x/9.0
+        colors[...,1] = y/9.0
+        colors[...,2] = z/9.0
+        ax.voxels(voxels, facecolors=colors)
+
+    @image_comparison(
+        baseline_images=['voxels-alpha'],
+        extensions=['png'],
+        remove_text=True
+    )
+    def test_alpha(self):
+        fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+        x, y, z = np.indices((10, 10, 10))
+        v1 = x == y
+        v2 = np.abs(x - y) < 2
+        voxels = v1 | v2
+        colors = np.zeros((10, 10, 10, 4))
+        colors[v2] = [1, 0, 0, 0.5]
+        colors[v1] = [0, 1, 0, 0.5]
+        v = ax.voxels(voxels, facecolors=colors)
+
+        assert type(v) is dict
+        for coord, poly in v.items():
+            assert voxels[coord], "faces returned for absent voxel"
+            assert isinstance(poly, art3d.Poly3DCollection)
+
+    @image_comparison(
+        baseline_images=['voxels-xyz'],
+        extensions=['png'],
+        tol=0.01
+    )
+    def test_xyz(self):
+        fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+        def midpoints(x):
+            sl = ()
+            for i in range(x.ndim):
+                x = (x[sl + np.index_exp[:-1]] +
+                     x[sl + np.index_exp[1:]]) / 2.0
+                sl += np.index_exp[:]
+            return x
+
+        # prepare some coordinates, and attach rgb values to each
+        r, g, b = np.indices((17, 17, 17)) / 16.0
+        rc = midpoints(r)
+        gc = midpoints(g)
+        bc = midpoints(b)
+
+        # define a sphere about [0.5, 0.5, 0.5]
+        sphere = (rc - 0.5)**2 + (gc - 0.5)**2 + (bc - 0.5)**2 < 0.5**2
+
+        # combine the color components
+        colors = np.zeros(sphere.shape + (3,))
+        colors[..., 0] = rc
+        colors[..., 1] = gc
+        colors[..., 2] = bc
+
+        # and plot everything
+        ax.voxels(r, g, b, sphere,
+                  facecolors=colors,
+                  edgecolors=np.clip(2*colors - 0.5, 0, 1),  # brighter
+                  linewidth=0.5)
+
+    def test_calling_conventions(self):
+        x, y, z = np.indices((3, 4, 5))
+        filled = np.ones((2, 3, 4))
+
+        fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+
+        # all the valid calling conventions
+        for kw in (dict(), dict(edgecolor='k')):
+            ax.voxels(filled, **kw)
+            ax.voxels(filled=filled, **kw)
+            ax.voxels(x, y, z, filled, **kw)
+            ax.voxels(x, y, z, filled=filled, **kw)
+
+        # duplicate argument
+        with pytest.raises(TypeError) as exc:
+            ax.voxels(x, y, z, filled, filled=filled)
+        exc.match(".*voxels.*")
+        # missing arguments
+        with pytest.raises(TypeError) as exc:
+            ax.voxels(x, y)
+        exc.match(".*voxels.*")
+        # x,y,z are positional only - this passes them on as attributes of
+        # Poly3DCollection
+        with pytest.raises(AttributeError):
+            ax.voxels(filled=filled, x=x, y=y, z=z)
+
+
 def test_inverted_cla():
     # Github PR #5450. Setting autoscale should reset
     # axes to be non-inverted.