diff --git a/doc/api/next_api_changes/behavior/31009-TK.rst b/doc/api/next_api_changes/behavior/31009-TK.rst
new file mode 100644
index 000000000000..1bd910f856ec
--- /dev/null
+++ b/doc/api/next_api_changes/behavior/31009-TK.rst
@@ -0,0 +1,6 @@
+imshow pixel alignment
+~~~~~~~~~~~~~~~~~~~~~~
+``imshow`` now rounds the output image dimensions to the nearest integer pixel
+rather than always using ``math.ceil``. This ensures that the edges of raster
+pixels align correctly with vector grid lines and patches, preventing 1-pixel
+overlaps or gaps when the image is small or pixels are large.
diff --git a/lib/matplotlib/image.py b/lib/matplotlib/image.py
index c1846f92608c..09505a910083 100644
--- a/lib/matplotlib/image.py
+++ b/lib/matplotlib/image.py
@@ -3,7 +3,6 @@
 operations.
 """
 
-import math
 import os
 import logging
 from pathlib import Path
@@ -420,16 +419,29 @@ def _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification=1.0,
                 .translate(-clipped_bbox.x0, -clipped_bbox.y0)
                 .scale(magnification)))
 
-        # So that the image is aligned with the edge of the Axes, we want to
-        # round up the output width to the next integer.  This also means
-        # scaling the transform slightly to account for the extra subpixel.
-        if ((not unsampled) and t.is_affine and round_to_pixel_border and
-                (out_width_base % 1.0 != 0.0 or out_height_base % 1.0 != 0.0)):
-            out_width = math.ceil(out_width_base)
-            out_height = math.ceil(out_height_base)
-            extra_width = (out_width - out_width_base) / out_width_base
-            extra_height = (out_height - out_height_base) / out_height_base
-            t += Affine2D().scale(1.0 + extra_width, 1.0 + extra_height)
+        # So that the image is aligned with the edge of the Axes, we round the
+        # output size to the nearest integer, and scale the transform slightly
+        # to account for any resulting subpixel difference.
+        if ((not unsampled) and t.is_affine and round_to_pixel_border):
+            out_width = round(out_width_base)
+            out_height = round(out_height_base)
+
+            if out_width > 0 and out_height > 0:
+                extra_width = (out_width - out_width_base) / out_width_base
+                extra_height = (out_height - out_height_base) / out_height_base
+                t += Affine2D().scale(1.0 + extra_width, 1.0 + extra_height)
+            else:
+                # Fallback for small positive base dimensions: ensure at least
+                # one pixel so we do not end up with a zero-sized output.
+                if out_width_base > 0:
+                    out_width = max(1, int(out_width_base))
+                else:
+                    out_width = int(out_width_base)
+
+                if out_height_base > 0:
+                    out_height = max(1, int(out_height_base))
+                else:
+                    out_height = int(out_height_base)
         else:
             out_width = int(out_width_base)
             out_height = int(out_height_base)
diff --git a/lib/matplotlib/tests/test_image.py b/lib/matplotlib/tests/test_image.py
index da7a198a2a94..b07ff7845ef1 100644
--- a/lib/matplotlib/tests/test_image.py
+++ b/lib/matplotlib/tests/test_image.py
@@ -1871,3 +1871,41 @@ def test_interpolation_stage_rgba_respects_alpha_param(fig_test, fig_ref, intp_s
             (im_rgb, new_array_alpha.reshape((ny, nx, 1))), axis=-1
         ), interpolation_stage=intp_stage
     )
+
+
+def test_imshow_pixel_rounding():
+    """
+    Test that imshow rounds output dimensions to the nearest integer
+    (matching grid snapping) rather than always ceiling them.
+    Regression test for: https://github.com/matplotlib/matplotlib/issues/31009
+    """
+    # 1. Setup a figure with known DPI
+    dpi = 100
+    fig = plt.figure(dpi=dpi)
+
+    # 2. Create an axes that occupies a precise fractional area
+    # We want the axes to be 10.4 pixels wide.
+    # width_inches = 10.4 / 100 = 0.104
+    fig.set_size_inches(1, 1)
+    ax = fig.add_axes([0, 0, 0.104, 0.104])  # [left, bottom, width, height]
+
+    # 3. Add an image that fills the axes
+    # Data is 1x1, Extent matches limits
+    ax.set_xlim(0, 1)
+    ax.set_ylim(0, 1)
+    im = ax.imshow([[1]], extent=[0, 1, 0, 1], interpolation='nearest')
+
+    # 4. Trigger the internal make_image call
+    # This invokes the logic we changed in image.py
+    fig.canvas.draw()
+    renderer = fig.canvas.get_renderer()
+
+    # im.make_image returns (image_array, x, y, transform)
+    # The image_array shape corresponds to the calculated pixel size
+    resampled_img, _, _, _ = im.make_image(renderer)
+    # 5. Assert the logic
+    # Theoretical width: 100 dpi * 1 inch * 0.104 fraction = 10.4 pixels.
+    # Old behavior (ceil): 11 pixels
+    # New behavior (round): 10 pixels
+    assert resampled_img.shape == (10, 10, 4), \
+        f"Expected 10x10 image (rounded), got {resampled_img.shape} (likely ceiled)"