diff --git a/Tests/test_image_resize.py b/Tests/test_image_resize.py
index d9ddf5009..1cd9c3800 100644
--- a/Tests/test_image_resize.py
+++ b/Tests/test_image_resize.py
@@ -44,9 +44,19 @@ class TestImagingCoreResize:
             self.resize(hopper("1"), (15, 12), Image.Resampling.BILINEAR)
         with pytest.raises(ValueError):
             self.resize(hopper("P"), (15, 12), Image.Resampling.BILINEAR)
-        with pytest.raises(ValueError):
-            self.resize(hopper("I;16"), (15, 12), Image.Resampling.BILINEAR)
-        for mode in ["L", "I", "F", "RGB", "RGBA", "CMYK", "YCbCr"]:
+        for mode in [
+            "L",
+            "I",
+            "I;16",
+            "I;16L",
+            "I;16B",
+            "I;16N",
+            "F",
+            "RGB",
+            "RGBA",
+            "CMYK",
+            "YCbCr",
+        ]:
             im = hopper(mode)
             r = self.resize(im, (15, 12), Image.Resampling.BILINEAR)
             assert r.mode == mode
diff --git a/src/_imaging.c b/src/_imaging.c
index 07d9a64cc..426036335 100644
--- a/src/_imaging.c
+++ b/src/_imaging.c
@@ -1579,16 +1579,12 @@ _putdata(ImagingObject *self, PyObject *args) {
                 int bigendian = 0;
                 if (image->type == IMAGING_TYPE_SPECIAL) {
                     // I;16*
-                    if (strcmp(image->mode, "I;16N") == 0) {
+                    if (strcmp(image->mode, "I;16B") == 0
 #ifdef WORDS_BIGENDIAN
-                        bigendian = 1;
-#else
-                        bigendian = 0;
+                        || strcmp(image->mode, "I;16N") == 0
 #endif
-                    } else if (strcmp(image->mode, "I;16B") == 0) {
+                    ) {
                         bigendian = 1;
-                    } else {
-                        bigendian = 0;
                     }
                 }
                 for (i = x = y = 0; i < n; i++) {
diff --git a/src/libImaging/Resample.c b/src/libImaging/Resample.c
index 222d6bca4..f5e386dc2 100644
--- a/src/libImaging/Resample.c
+++ b/src/libImaging/Resample.c
@@ -460,6 +460,83 @@ ImagingResampleVertical_8bpc(
     ImagingSectionLeave(&cookie);
 }
 
+void
+ImagingResampleHorizontal_16bpc(
+    Imaging imOut, Imaging imIn, int offset, int ksize, int *bounds, double *kk
+) {
+    ImagingSectionCookie cookie;
+    double ss;
+    int xx, yy, x, xmin, xmax, ss_int;
+    double *k;
+
+    int bigendian = 0;
+    if (strcmp(imIn->mode, "I;16N") == 0
+#ifdef WORDS_BIGENDIAN
+        || strcmp(imIn->mode, "I;16B") == 0
+#endif
+    ) {
+        bigendian = 1;
+    }
+
+    ImagingSectionEnter(&cookie);
+    for (yy = 0; yy < imOut->ysize; yy++) {
+        for (xx = 0; xx < imOut->xsize; xx++) {
+            xmin = bounds[xx * 2 + 0];
+            xmax = bounds[xx * 2 + 1];
+            k = &kk[xx * ksize];
+            ss = 0.0;
+            for (x = 0; x < xmax; x++) {
+                ss += (imIn->image8[yy + offset][(x + xmin) * 2 + (bigendian ? 1 : 0)] +
+                       (imIn->image8[yy + offset][(x + xmin) * 2 + (bigendian ? 0 : 1)]
+                        << 8)) *
+                      k[x];
+            }
+            ss_int = ROUND_UP(ss);
+            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = CLIP8(ss_int % 256);
+            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = CLIP8(ss_int >> 8);
+        }
+    }
+    ImagingSectionLeave(&cookie);
+}
+
+void
+ImagingResampleVertical_16bpc(
+    Imaging imOut, Imaging imIn, int offset, int ksize, int *bounds, double *kk
+) {
+    ImagingSectionCookie cookie;
+    double ss;
+    int xx, yy, y, ymin, ymax, ss_int;
+    double *k;
+
+    int bigendian = 0;
+    if (strcmp(imIn->mode, "I;16N") == 0
+#ifdef WORDS_BIGENDIAN
+        || strcmp(imIn->mode, "I;16B") == 0
+#endif
+    ) {
+        bigendian = 1;
+    }
+
+    ImagingSectionEnter(&cookie);
+    for (yy = 0; yy < imOut->ysize; yy++) {
+        ymin = bounds[yy * 2 + 0];
+        ymax = bounds[yy * 2 + 1];
+        k = &kk[yy * ksize];
+        for (xx = 0; xx < imOut->xsize; xx++) {
+            ss = 0.0;
+            for (y = 0; y < ymax; y++) {
+                ss += (imIn->image8[y + ymin][xx * 2 + (bigendian ? 1 : 0)] +
+                       (imIn->image8[y + ymin][xx * 2 + (bigendian ? 0 : 1)] << 8)) *
+                      k[y];
+            }
+            ss_int = ROUND_UP(ss);
+            imOut->image8[yy][xx * 2 + (bigendian ? 1 : 0)] = CLIP8(ss_int % 256);
+            imOut->image8[yy][xx * 2 + (bigendian ? 0 : 1)] = CLIP8(ss_int >> 8);
+        }
+    }
+    ImagingSectionLeave(&cookie);
+}
+
 void
 ImagingResampleHorizontal_32bpc(
     Imaging imOut, Imaging imIn, int offset, int ksize, int *bounds, double *kk
@@ -574,7 +651,12 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter, float box[4]) {
     }
 
     if (imIn->type == IMAGING_TYPE_SPECIAL) {
-        return (Imaging)ImagingError_ModeError();
+        if (strncmp(imIn->mode, "I;16", 4) == 0) {
+            ResampleHorizontal = ImagingResampleHorizontal_16bpc;
+            ResampleVertical = ImagingResampleVertical_16bpc;
+        } else {
+            return (Imaging)ImagingError_ModeError();
+        }
     } else if (imIn->image8) {
         ResampleHorizontal = ImagingResampleHorizontal_8bpc;
         ResampleVertical = ImagingResampleVertical_8bpc;