diff --git a/lib/mpl_toolkits/mplot3d/axes3d.py b/lib/mpl_toolkits/mplot3d/axes3d.py
@@ -1437,17 +1437,6 @@ def plot_surface(self, X, Y, Z, *args, norm=None, vmin=None,
            the input data is larger, it will be downsampled (by slicing) to
            these numbers of points.
 
-        .. note::
-
-           To maximize rendering speed consider setting *rstride* and *cstride*
-           to divisors of the number of rows minus 1 and columns minus 1
-           respectively. For example, given 51 rows rstride can be any of the
-           divisors of 50.
-
-           Similarly, a setting of *rstride* and *cstride* equal to 1 (or
-           *rcount* and *ccount* equal the number of rows and columns) can use
-           the optimized path.
-
         Parameters
         ----------
         X, Y, Z : 2d arrays
@@ -1533,7 +1522,7 @@ def plot_surface(self, X, Y, Z, *args, norm=None, vmin=None,
             cstride = int(max(np.ceil(cols / ccount), 1))
 
         if 'facecolors' in kwargs:
-            fcolors = kwargs.pop('facecolors')
+            fcolors = np.asarray(kwargs.pop('facecolors'))
         else:
             color = kwargs.pop('color', None)
             if color is None:
@@ -1551,33 +1540,31 @@ def plot_surface(self, X, Y, Z, *args, norm=None, vmin=None,
                         "semantic or raise an error in matplotlib 3.3. "
                         "Please use shade=False instead.")
 
-        colset = []  # the sampled facecolor
-        if (rows - 1) % rstride == 0 and \
-           (cols - 1) % cstride == 0 and \
-           fcolors is None:
-            polys = np.stack(
-                [cbook._array_patch_perimeters(a, rstride, cstride)
-                 for a in (X, Y, Z)],
-                axis=-1)
-        else:
-            # evenly spaced, and including both endpoints
-            row_inds = list(range(0, rows-1, rstride)) + [rows-1]
-            col_inds = list(range(0, cols-1, cstride)) + [cols-1]
-
-            polys = []
-            for rs, rs_next in zip(row_inds[:-1], row_inds[1:]):
-                for cs, cs_next in zip(col_inds[:-1], col_inds[1:]):
-                    ps = [
-                        # +1 ensures we share edges between polygons
-                        cbook._array_perimeter(a[rs:rs_next+1, cs:cs_next+1])
-                        for a in (X, Y, Z)
-                    ]
-                    # ps = np.stack(ps, axis=-1)
-                    ps = np.array(ps).T
-                    polys.append(ps)
-
-                    if fcolors is not None:
-                        colset.append(fcolors[rs][cs])
+        # Calculate the minimal amount of padding that will allow us to extract
+        # vectorized patches.
+        rrem = (rows - 1) % rstride
+        rpadding = rstride - rrem if rrem != 0 else 0
+        crem = (cols - 1) % cstride
+        cpadding = cstride - crem if crem != 0 else 0
+
+        def pad_for_patches(a):
+            if cpadding == 0 and rpadding == 0:
+                return a
+            result = np.empty_like(a, shape=(rows + rpadding, cols + cpadding))
+            result[:rows, :cols] = a
+            result[:rows, cols:] = a[:, -1:]
+            result[rows:, :cols] = a[-1:, :]
+            result[rows:, cols:] = a[-1, -1]
+            return result
+
+        polys = np.stack(
+            [cbook._array_patch_perimeters(pad_for_patches(a),
+                                           rstride, cstride)
+                for a in (X, Y, Z)],
+            axis=-1)
+        if fcolors is not None:
+            colset = fcolors[np.ix_(np.arange(0, rows-1, rstride),
+                                    np.arange(0, cols-1, cstride))]
 
         # note that the striding causes some polygons to have more coordinates
         # than others

diff --git a/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/surface3d_shaded.png b/lib/mpl_toolkits/tests/baseline_images/test_mplot3d/surface3d_shaded.png