diff --git a/Tests/test_image_resample.py b/Tests/test_image_resample.py
index 73b25ed51..d4bd81559 100644
--- a/Tests/test_image_resample.py
+++ b/Tests/test_image_resample.py
@@ -165,6 +165,32 @@ class TestImagingCoreResampleAccuracy:
         for channel in case.split():
             self.check_case(channel, self.make_sample(data, (8, 8)))
 
+    @pytest.mark.parametrize("mode", ("RGBX", "RGB", "La", "L"))
+    def test_reduce_mks2013(self, mode: str) -> None:
+        case = self.make_case(mode, (16, 16), 0xE1)
+        case = case.resize((8, 8), Image.Resampling.MKS2013)
+        # fmt: off
+        data = ("e1 e1 e9 dc"
+                "e1 e1 e9 dc"
+                "e9 e9 f1 e3"
+                "dc dc e4 d8")
+        # fmt: on
+        for channel in case.split():
+            self.check_case(channel, self.make_sample(data, (8, 8)))
+
+    @pytest.mark.parametrize("mode", ("RGBX", "RGB", "La", "L"))
+    def test_reduce_mks2021(self, mode: str) -> None:
+        case = self.make_case(mode, (16, 16), 0xE1)
+        case = case.resize((8, 8), Image.Resampling.MKS2021)
+        # fmt: off
+        data = ("e1 e1 e3 d7"
+                "e1 e1 e3 d7"
+                "e3 e3 e5 d9"
+                "d7 d7 d9 ce")
+        # fmt: on
+        for channel in case.split():
+            self.check_case(channel, self.make_sample(data, (8, 8)))
+
     @pytest.mark.parametrize("mode", ("RGBX", "RGB", "La", "L"))
     def test_enlarge_box(self, mode: str) -> None:
         case = self.make_case(mode, (2, 2), 0xE1)
@@ -226,6 +252,36 @@ class TestImagingCoreResampleAccuracy:
         for channel in case.split():
             self.check_case(channel, self.make_sample(data, (12, 12)))
 
+    @pytest.mark.parametrize("mode", ("RGBX", "RGB", "La", "L"))
+    def test_enlarge_mks2013(self, mode: str) -> None:
+        case = self.make_case(mode, (6, 6), 0xE1)
+        case = case.resize((12, 12), Image.Resampling.MKS2013)
+        data = (
+            "e1 e1 e2 ef fb be"
+            "e1 e1 e2 ef fb be"
+            "e2 e2 e3 f1 fd bf"
+            "ef ef f0 ff ff c7"
+            "fb fb fc ff ff cf"
+            "be be bf c7 cf a8"
+        )
+        for channel in case.split():
+            self.check_case(channel, self.make_sample(data, (12, 12)))
+
+    @pytest.mark.parametrize("mode", ("RGBX", "RGB", "La", "L"))
+    def test_enlarge_mks2021(self, mode: str) -> None:
+        case = self.make_case(mode, (6, 6), 0xE1)
+        case = case.resize((12, 12), Image.Resampling.MKS2021)
+        data = (
+            "e3 e1 df e9 f5 bb"
+            "e1 df dd e7 f3 b9"
+            "df dd db e5 f1 b8"
+            "e9 e7 e5 ef fc be"
+            "f5 f3 f0 fc ff c5"
+            "bb ba b8 bf c6 a3"
+        )
+        for channel in case.split():
+            self.check_case(channel, self.make_sample(data, (12, 12)))
+
     def test_box_filter_correct_range(self) -> None:
         im = Image.new("RGB", (8, 8), "#1688ff").resize(
             (100, 100), Image.Resampling.BOX
@@ -309,6 +365,8 @@ class TestCoreResampleAlphaCorrect:
         self.run_levels_case(case.resize((512, 32), Image.Resampling.HAMMING))
         self.run_levels_case(case.resize((512, 32), Image.Resampling.BICUBIC))
         self.run_levels_case(case.resize((512, 32), Image.Resampling.LANCZOS))
+        self.run_levels_case(case.resize((512, 32), Image.Resampling.MKS2013))
+        self.run_levels_case(case.resize((512, 32), Image.Resampling.MKS2021))
 
     @pytest.mark.xfail(reason="Current implementation isn't precise enough")
     def test_levels_la(self) -> None:
@@ -318,6 +376,8 @@ class TestCoreResampleAlphaCorrect:
         self.run_levels_case(case.resize((512, 32), Image.Resampling.HAMMING))
         self.run_levels_case(case.resize((512, 32), Image.Resampling.BICUBIC))
         self.run_levels_case(case.resize((512, 32), Image.Resampling.LANCZOS))
+        self.run_levels_case(case.resize((512, 32), Image.Resampling.MKS2013))
+        self.run_levels_case(case.resize((512, 32), Image.Resampling.MKS2021))
 
     def make_dirty_case(
         self, mode: str, clean_pixel: tuple[int, ...], dirty_pixel: tuple[int, ...]
@@ -360,6 +420,12 @@ class TestCoreResampleAlphaCorrect:
         self.run_dirty_case(
             case.resize((20, 20), Image.Resampling.LANCZOS), (255, 255, 0)
         )
+        self.run_dirty_case(
+            case.resize((20, 20), Image.Resampling.MKS2013), (255, 255, 0)
+        )
+        self.run_dirty_case(
+            case.resize((20, 20), Image.Resampling.MKS2021), (255, 255, 0)
+        )
 
     def test_dirty_pixels_la(self) -> None:
         case = self.make_dirty_case("LA", (255, 128), (0, 0))
@@ -368,6 +434,8 @@ class TestCoreResampleAlphaCorrect:
         self.run_dirty_case(case.resize((20, 20), Image.Resampling.HAMMING), (255,))
         self.run_dirty_case(case.resize((20, 20), Image.Resampling.BICUBIC), (255,))
         self.run_dirty_case(case.resize((20, 20), Image.Resampling.LANCZOS), (255,))
+        self.run_dirty_case(case.resize((20, 20), Image.Resampling.MKS2013), (255,))
+        self.run_dirty_case(case.resize((20, 20), Image.Resampling.MKS2021), (255,))
 
 
 class TestCoreResamplePasses:
@@ -453,6 +521,8 @@ class TestCoreResampleBox:
             Image.Resampling.HAMMING,
             Image.Resampling.BICUBIC,
             Image.Resampling.LANCZOS,
+            Image.Resampling.MKS2013,
+            Image.Resampling.MKS2021,
         ),
     )
     def test_wrong_arguments(self, resample: Image.Resampling) -> None:
diff --git a/Tests/test_image_resize.py b/Tests/test_image_resize.py
index 270500a44..9e549cc71 100644
--- a/Tests/test_image_resize.py
+++ b/Tests/test_image_resize.py
@@ -72,6 +72,8 @@ class TestImagingCoreResize:
             Image.Resampling.HAMMING,
             Image.Resampling.BICUBIC,
             Image.Resampling.LANCZOS,
+            Image.Resampling.MKS2013,
+            Image.Resampling.MKS2021,
         ),
     )
     def test_reduce_filters(self, resample: Image.Resampling) -> None:
@@ -88,6 +90,8 @@ class TestImagingCoreResize:
             Image.Resampling.HAMMING,
             Image.Resampling.BICUBIC,
             Image.Resampling.LANCZOS,
+            Image.Resampling.MKS2013,
+            Image.Resampling.MKS2021,
         ),
     )
     def test_enlarge_filters(self, resample: Image.Resampling) -> None:
@@ -104,6 +108,8 @@ class TestImagingCoreResize:
             Image.Resampling.HAMMING,
             Image.Resampling.BICUBIC,
             Image.Resampling.LANCZOS,
+            Image.Resampling.MKS2013,
+            Image.Resampling.MKS2021,
         ),
     )
     @pytest.mark.parametrize(
@@ -154,6 +160,8 @@ class TestImagingCoreResize:
             Image.Resampling.HAMMING,
             Image.Resampling.BICUBIC,
             Image.Resampling.LANCZOS,
+            Image.Resampling.MKS2013,
+            Image.Resampling.MKS2021,
         ),
     )
     def test_enlarge_zero(self, resample: Image.Resampling) -> None:
diff --git a/docs/handbook/concepts.rst b/docs/handbook/concepts.rst
index 46f612be3..0a4801af6 100644
--- a/docs/handbook/concepts.rst
+++ b/docs/handbook/concepts.rst
@@ -217,6 +217,23 @@ pixel, the Python Imaging Library provides different resampling *filters*.
 
     .. versionadded:: 1.1.3
 
+.. data:: Resampling.MKS2013
+    :noindex:
+
+    Calculate the output pixel value using the Magic Kernel Sharp 2013 filter
+    (a quadratic B-spline composed with a sharpening kernel) on all pixels that
+    may contribute to the output value. This filter can only be used with the
+    :py:meth:`~PIL.Image.Image.resize` and :py:meth:`~PIL.Image.Image.thumbnail`
+    methods.
+
+.. data:: Resampling.MKS2021
+    :noindex:
+
+    Calculate the output pixel value using the Magic Kernel Sharp 2021 filter
+    (a quadratic B-spline composed with a sharpening kernel) on all pixels that
+    may contribute to the output value. This filter can only be used with the
+    :py:meth:`~PIL.Image.Image.resize` and :py:meth:`~PIL.Image.Image.thumbnail`
+    methods.
 
 Filters comparison table
 ~~~~~~~~~~~~~~~~~~~~~~~~
@@ -237,3 +254,7 @@ Filters comparison table
 +---------------------------+-------------+-----------+-------------+
 |:data:`Resampling.LANCZOS` | ⭐⭐⭐⭐    | ⭐⭐⭐⭐  | ⭐          |
 +---------------------------+-------------+-----------+-------------+
+|:data:`Resampling.MKS2013` | ⭐⭐⭐⭐    | ⭐⭐⭐⭐  | ⭐          |
++---------------------------+-------------+-----------+-------------+
+|:data:`Resampling.MKS2021` | ⭐⭐⭐⭐⭐  | ⭐⭐⭐⭐⭐| ⭐          |
++---------------------------+-------------+-----------+-------------+
diff --git a/src/PIL/Image.py b/src/PIL/Image.py
index 262b5478b..e1610272f 100644
--- a/src/PIL/Image.py
+++ b/src/PIL/Image.py
@@ -147,6 +147,8 @@ class Resampling(IntEnum):
     HAMMING = 5
     BICUBIC = 3
     LANCZOS = 1
+    MKS2013 = 6
+    MKS2021 = 7
 
 
 _filters_support = {
@@ -155,6 +157,8 @@ _filters_support = {
     Resampling.HAMMING: 1.0,
     Resampling.BICUBIC: 2.0,
     Resampling.LANCZOS: 3.0,
+    Resampling.MKS2013: 2.5,
+    Resampling.MKS2021: 4.5,
 }
 
 
@@ -2211,7 +2215,8 @@ class Image:
         :param resample: An optional resampling filter.  This can be
            one of :py:data:`Resampling.NEAREST`, :py:data:`Resampling.BOX`,
            :py:data:`Resampling.BILINEAR`, :py:data:`Resampling.HAMMING`,
-           :py:data:`Resampling.BICUBIC` or :py:data:`Resampling.LANCZOS`.
+           :py:data:`Resampling.BICUBIC`, :py:data:`Resampling.LANCZOS`,
+           :py:data:`Resampling.MKS2013`, or :py:data:`Resampling.MKS2021`.
            If the image has mode "1" or "P", it is always set to
            :py:data:`Resampling.NEAREST`. Otherwise, the default filter is
            :py:data:`Resampling.BICUBIC`. See: :ref:`concept-filters`.
@@ -2243,6 +2248,8 @@ class Image:
             Resampling.LANCZOS,
             Resampling.BOX,
             Resampling.HAMMING,
+            Resampling.MKS2013,
+            Resampling.MKS2021,
         ):
             msg = f"Unknown resampling filter ({resample})."
 
@@ -2255,6 +2262,8 @@ class Image:
                     (Resampling.BICUBIC, "Image.Resampling.BICUBIC"),
                     (Resampling.BOX, "Image.Resampling.BOX"),
                     (Resampling.HAMMING, "Image.Resampling.HAMMING"),
+                    (Resampling.MKS2013, "Image.Resampling.MKS2013"),
+                    (Resampling.MKS2021, "Image.Resampling.MKS2021"),
                 )
             ]
             msg += f" Use {', '.join(filters[:-1])} or {filters[-1]}"
@@ -2706,7 +2715,8 @@ class Image:
         :param resample: Optional resampling filter.  This can be one
            of :py:data:`Resampling.NEAREST`, :py:data:`Resampling.BOX`,
            :py:data:`Resampling.BILINEAR`, :py:data:`Resampling.HAMMING`,
-           :py:data:`Resampling.BICUBIC` or :py:data:`Resampling.LANCZOS`.
+           :py:data:`Resampling.BICUBIC`, :py:data:`Resampling.LANCZOS`,
+           :py:data:`Resampling.MKS2013`, or :py:data:`Resampling.MKS2021`.
            If omitted, it defaults to :py:data:`Resampling.BICUBIC`.
            (was :py:data:`Resampling.NEAREST` prior to version 2.5.0).
            See: :ref:`concept-filters`.
@@ -2920,11 +2930,19 @@ class Image:
             Resampling.BILINEAR,
             Resampling.BICUBIC,
         ):
-            if resample in (Resampling.BOX, Resampling.HAMMING, Resampling.LANCZOS):
+            if resample in (
+                Resampling.BOX,
+                Resampling.HAMMING,
+                Resampling.LANCZOS,
+                Resampling.MKS2013,
+                Resampling.MKS2021,
+            ):
                 unusable: dict[int, str] = {
                     Resampling.BOX: "Image.Resampling.BOX",
                     Resampling.HAMMING: "Image.Resampling.HAMMING",
                     Resampling.LANCZOS: "Image.Resampling.LANCZOS",
+                    Resampling.MKS2013: "Image.Resampling.MKS2013",
+                    Resampling.MKS2021: "Image.Resampling.MKS2021",
                 }
                 msg = unusable[resample] + f" ({resample}) cannot be used."
             else:
diff --git a/src/libImaging/Imaging.h b/src/libImaging/Imaging.h
index bfe67d462..675a68046 100644
--- a/src/libImaging/Imaging.h
+++ b/src/libImaging/Imaging.h
@@ -294,6 +294,8 @@ ImagingError_ValueError(const char *message);
 #define IMAGING_TRANSFORM_HAMMING 5
 #define IMAGING_TRANSFORM_BICUBIC 3
 #define IMAGING_TRANSFORM_LANCZOS 1
+#define IMAGING_TRANSFORM_MKS2013 6
+#define IMAGING_TRANSFORM_MKS2021 7
 
 typedef int (*ImagingTransformMap)(double *X, double *Y, int x, int y, void *data);
 typedef int (*ImagingTransformFilter)(void *out, Imaging im, double x, double y);
diff --git a/src/libImaging/Resample.c b/src/libImaging/Resample.c
index b114e0023..ae0186c71 100644
--- a/src/libImaging/Resample.c
+++ b/src/libImaging/Resample.c
@@ -79,11 +79,55 @@ lanczos_filter(double x) {
     return 0.0;
 }
 
+static inline double
+mks_2013_filter(double x) {
+    /* https://johncostella.com/magic/ */
+    if (x < 0) {
+        x = -x;
+    }
+    if (x < 0.5) {
+        return 17.0 / 16.0 - 7.0 / 4.0 * pow(x, 2);
+    }
+    if (x < 1.5) {
+        return (1.0 - x) * (7.0 / 4.0 - x);
+    }
+    if (x < 2.5) {
+        return -1.0 / 8.0 * pow(x - 5.0 / 2.0, 2);
+    }
+    return 0;
+}
+
+static inline double
+mks_2021_filter(double x) {
+    /* https://johncostella.com/magic/ */
+    if (x < 0) {
+        x = -x;
+    }
+    if (x < 0.5) {
+        return 577.0 / 576.0 - 239.0 / 144.0 * pow(x, 2);
+    }
+    if (x < 1.5) {
+        return 35.0 / 36.0 * (x - 1.0) * (x - 239.0 / 140.0);
+    }
+    if (x < 2.5) {
+        return 1.0 / 6.0 * (x - 2.0) * (65.0 / 24.0 - x);
+    }
+    if (x < 3.5) {
+        return 1.0 / 36.0 * (x - 3.0) * (x - 15.0 / 4.0);
+    }
+    if (x < 4.5) {
+        return -1.0 / 288.0 * pow(x - 9.0 / 2.0, 2);
+    }
+    return 0;
+}
+
 static struct filter BOX = {box_filter, 0.5};
 static struct filter BILINEAR = {bilinear_filter, 1.0};
 static struct filter HAMMING = {hamming_filter, 1.0};
 static struct filter BICUBIC = {bicubic_filter, 2.0};
 static struct filter LANCZOS = {lanczos_filter, 3.0};
+static struct filter MKS2013 = {mks_2013_filter, 2.5};
+static struct filter MKS2021 = {mks_2021_filter, 4.5};
 
 /* 8 bits for result. Filter can have negative areas.
    In one cases the sum of the coefficients will be negative,
@@ -695,6 +739,12 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter, float box[4]) {
         case IMAGING_TRANSFORM_LANCZOS:
             filterp = &LANCZOS;
             break;
+        case IMAGING_TRANSFORM_MKS2013:
+            filterp = &MKS2013;
+            break;
+        case IMAGING_TRANSFORM_MKS2021:
+            filterp = &MKS2021;
+            break;
         default:
             return (Imaging)ImagingError_ValueError("unsupported resampling filter");
     }