Make .quantize(palette=...) use exact colors when possible, and improve dithering quality

Tests/benchmarks.py

@@ -898,3 +898,5 @@ def test_quantize_to_palette(
         result = bench(lambda: im._new(im.im.convert(output_mode, dither, palette.im)))
     assert result.mode == output_mode
     benchmark_save(result)
+    if palette_type == "exact":
+        assert result.convert("RGB").tobytes() == im.tobytes()

Tests/test_image_entropy.py

@@ -13,7 +13,7 @@ def entropy(mode: str) -> float:
     assert entropy("L") == pytest.approx(7.063008716585465)
     assert entropy("I") == pytest.approx(7.063008716585465)
     assert entropy("F") == pytest.approx(7.063008716585465)
-    assert entropy("P") == pytest.approx(5.082506854662517)
+    assert entropy("P") == pytest.approx(5.085630558679487)
     assert entropy("RGB") == pytest.approx(8.821286587714319)
     assert entropy("RGBA") == pytest.approx(7.42724306524488)
     assert entropy("CMYK") == pytest.approx(7.4272430652448795)

Tests/test_image_histogram.py

@@ -14,7 +14,7 @@ def histogram(mode: str) -> tuple[int, int, int]:
     assert histogram("La") == (512, 0, 16384)
     assert histogram("I") == (256, 0, 662)
     assert histogram("F") == (256, 0, 662)
-    assert histogram("P") == (256, 0, 1551)
+    assert histogram("P") == (256, 0, 1557)
     assert histogram("PA") == (512, 0, 16384)
     assert histogram("RGB") == (768, 4, 675)
     assert histogram("RGBA") == (1024, 0, 16384)

Tests/test_image_quantize.py

@@ -94,6 +94,43 @@ def test_quantize_no_dither2() -> None:
         assert px[x, 0] == (0 if x < 5 else 1)
 
 
+@pytest.mark.parametrize("dither", (Image.Dither.NONE, Image.Dither.FLOYDSTEINBERG))
+def test_quantize_exact_palette_matches(dither: Image.Dither) -> None:
+    # These colors share a slot in the nearest-color cache.
+    colors = ((252, 252, 252), (255, 255, 255))
+    im = Image.new("RGB", (2, 1))
+    im.putdata(colors)
+
+    palette = Image.new("P", (1, 1))
+    palette.putpalette(tuple(channel for color in colors for channel in color))
+    quantized = im.quantize(palette=palette, dither=dither)
+
+    assert quantized.tobytes() == b"\x00\x01"
+    assert quantized.convert("RGB").tobytes() == im.tobytes()
+
+
+def test_quantize_exact_palette_matches_maximum_colors() -> None:
+    colors = tuple((r, 0, 0) for r in range(256))
+    im = Image.new("RGB", (256, 1))
+    im.putdata(colors)
+
+    palette = Image.new("P", (1, 1))
+    palette.putpalette(tuple(channel for color in colors for channel in color))
+    quantized = im.quantize(palette=palette, dither=Image.Dither.NONE)
+
+    assert quantized.tobytes() == bytes(range(256))
+    assert quantized.convert("RGB").tobytes() == im.tobytes()
+
+
+def test_quantize_duplicate_exact_palette_color_uses_first_index() -> None:
+    color = (17, 34, 51)
+    palette = Image.new("P", (1, 1))
+    palette.putpalette(color * 2)
+
+    im = Image.new("RGB", (1, 1), color)
+    assert im.quantize(palette=palette, dither=Image.Dither.NONE).getpixel((0, 0)) == 0
+
+
 def test_quantize_dither_diff() -> None:
     image = hopper()
     with Image.open("Tests/images/caption_6_33_22.png") as palette:

setup.py

@@ -88,6 +88,7 @@ def get_version() -> str:
     "Chops",
     "ColorLUT",
     "Convert",
+    "ConvertToPalette",
     "ConvertYCbCr",
     "Copy",
     "Crop",

src/PIL/Image.py

@@ -1298,8 +1298,9 @@ def quantize(
                        and :data:`Quantize.MAXCOVERAGE` do not support RGBA images, so
                        :data:`Quantize.FASTOCTREE` is used by default instead.
         :param kmeans: Integer greater than or equal to zero.
-        :param palette: Quantize to the palette of given
-                        :py:class:`PIL.Image.Image`.
+        :param palette: Quantize to the palette of given :py:class:`PIL.Image.Image`.
+                        The `colors`, `method` and `kmeans` parameters are ignored
+                        if a reference palette is used.
         :param dither: Dithering method, used when converting from
            mode "RGB" to "P" or from "RGB" or "L" to "1".
            Available methods are :data:`Dither.NONE` or :data:`Dither.FLOYDSTEINBERG`

src/_imaging.c

@@ -365,6 +365,10 @@ static const char *no_palette = "image has no palette";
 static const char *readonly = "image is readonly";
 /* static const char* no_content = "image has no content"; */
 
+/**
+ * Set a MemoryError exception and return NULL.
+ * @return Always NULL.
+ */
 void *
 ImagingError_MemoryError(void) {
     return PyErr_NoMemory();

src/libImaging/Convert.c

@@ -33,6 +33,7 @@
  */
 
 #include "Imaging.h"
+#include "ConvertToPalette.h"
 
 /* ITU-R Recommendation 601-2 (assuming nonlinear RGB) */
 #define L(rgb) ((INT32)(rgb)[0] * 299 + (INT32)(rgb)[1] * 587 + (INT32)(rgb)[2] * 114)
@@ -1141,206 +1142,6 @@ frompalette(Imaging imOut, Imaging imIn, const ModeID mode) {
 #if defined(_MSC_VER)
 #pragma optimize("", off)
 #endif
-static Imaging
-topalette(
-    Imaging imOut, Imaging imIn, const ModeID mode, ImagingPalette inpalette, int dither
-) {
-    ImagingSectionCookie cookie;
-    int alpha;
-    int x, y;
-    ImagingPalette palette = inpalette;
-
-    /* Map L or RGB/RGBX/RGBA/RGBa to palette image */
-    if (imIn->mode != IMAGING_MODE_L && imIn->mode != IMAGING_MODE_RGB &&
-        imIn->mode != IMAGING_MODE_RGBX && imIn->mode != IMAGING_MODE_RGBA &&
-        imIn->mode != IMAGING_MODE_RGBa) {
-        return (Imaging)ImagingError_ValueError("conversion not supported");
-    }
-
-    alpha = mode == IMAGING_MODE_PA;
-
-    if (palette == NULL) {
-        /* FIXME: make user configurable */
-        if (imIn->bands == 1) {
-            palette = ImagingPaletteNew(IMAGING_MODE_RGB);
-
-            palette->size = 256;
-            int i;
-            for (i = 0; i < 256; i++) {
-                palette->palette[i * 4] = palette->palette[i * 4 + 1] =
-                    palette->palette[i * 4 + 2] = (UINT8)i;
-            }
-        } else {
-            palette = ImagingPaletteNewBrowser(); /* Standard colour cube */
-        }
-    }
-
-    if (!palette) {
-        return (Imaging)ImagingError_ValueError("no palette");
-    }
-
-    imOut = ImagingNew2Dirty(mode, imOut, imIn);
-    if (!imOut) {
-        if (palette != inpalette) {
-            ImagingPaletteDelete(palette);
-        }
-        return NULL;
-    }
-
-    ImagingPaletteDelete(imOut->palette);
-    imOut->palette = ImagingPaletteDuplicate(palette);
-
-    if (imIn->bands == 1) {
-        /* grayscale image */
-
-        /* Grayscale palette: copy data as is */
-        ImagingSectionEnter(&cookie);
-        for (y = 0; y < imIn->ysize; y++) {
-            if (alpha) {
-                l2rgb((UINT8 *)imOut->image[y], (UINT8 *)imIn->image[y], imIn->xsize);
-            } else {
-                memcpy(imOut->image[y], imIn->image[y], imIn->linesize);
-            }
-        }
-        ImagingSectionLeave(&cookie);
-
-    } else {
-        /* colour image */
-
-        /* Create mapping cache */
-        if (ImagingPaletteCachePrepare(palette) < 0) {
-            ImagingDelete(imOut);
-            if (palette != inpalette) {
-                ImagingPaletteDelete(palette);
-            }
-            return NULL;
-        }
-
-        if (dither) {
-            /* floyd-steinberg dither */
-
-            int *errors;
-            errors = calloc(imIn->xsize + 1, sizeof(int) * 3);
-            if (!errors) {
-                ImagingDelete(imOut);
-                return ImagingError_MemoryError();
-            }
-
-            /* Map each pixel to the nearest palette entry */
-            ImagingSectionEnter(&cookie);
-            for (y = 0; y < imIn->ysize; y++) {
-                int r, r0, r1, r2;
-                int g, g0, g1, g2;
-                int b, b0, b1, b2;
-                UINT8 *in = (UINT8 *)imIn->image[y];
-                UINT8 *out = alpha ? (UINT8 *)imOut->image32[y] : imOut->image8[y];
-                int *e = errors;
-
-                r = r0 = r1 = 0;
-                g = g0 = g1 = 0;
-                b = b0 = b1 = b2 = 0;
-
-                for (x = 0; x < imIn->xsize; x++, in += 4) {
-                    int d2;
-                    INT16 *cache;
-
-                    r = CLIP8(in[0] + (r + e[3 + 0]) / 16);
-                    g = CLIP8(in[1] + (g + e[3 + 1]) / 16);
-                    b = CLIP8(in[2] + (b + e[3 + 2]) / 16);
-
-                    /* get closest colour */
-                    cache = &ImagingPaletteCache(palette, r, g, b);
-                    if (cache[0] == 0x100) {
-                        ImagingPaletteCacheUpdate(palette, r, g, b);
-                    }
-                    if (alpha) {
-                        out[x * 4] = out[x * 4 + 1] = out[x * 4 + 2] = (UINT8)cache[0];
-                        out[x * 4 + 3] = 255;
-                    } else {
-                        out[x] = (UINT8)cache[0];
-                    }
-
-                    r -= (int)palette->palette[cache[0] * 4];
-                    g -= (int)palette->palette[cache[0] * 4 + 1];
-                    b -= (int)palette->palette[cache[0] * 4 + 2];
-
-                    /* propagate errors (don't ask ;-) */
-                    r2 = r;
-                    d2 = r + r;
-                    r += d2;
-                    e[0] = r + r0;
-                    r += d2;
-                    r0 = r + r1;
-                    r1 = r2;
-                    r += d2;
-                    g2 = g;
-                    d2 = g + g;
-                    g += d2;
-                    e[1] = g + g0;
-                    g += d2;
-                    g0 = g + g1;
-                    g1 = g2;
-                    g += d2;
-                    b2 = b;
-                    d2 = b + b;
-                    b += d2;
-                    e[2] = b + b0;
-                    b += d2;
-                    b0 = b + b1;
-                    b1 = b2;
-                    b += d2;
-
-                    e += 3;
-                }
-
-                e[0] = b0;
-                e[1] = b1;
-                e[2] = b2;
-            }
-            ImagingSectionLeave(&cookie);
-            free(errors);
-
-        } else {
-            /* closest colour */
-            ImagingSectionEnter(&cookie);
-            for (y = 0; y < imIn->ysize; y++) {
-                int r, g, b;
-                UINT8 *in = (UINT8 *)imIn->image[y];
-                UINT8 *out = alpha ? (UINT8 *)imOut->image32[y] : imOut->image8[y];
-
-                for (x = 0; x < imIn->xsize; x++, in += 4) {
-                    INT16 *cache;
-
-                    r = in[0];
-                    g = in[1];
-                    b = in[2];
-
-                    /* get closest colour */
-                    cache = &ImagingPaletteCache(palette, r, g, b);
-                    if (cache[0] == 0x100) {
-                        ImagingPaletteCacheUpdate(palette, r, g, b);
-                    }
-                    if (alpha) {
-                        out[x * 4] = out[x * 4 + 1] = out[x * 4 + 2] = (UINT8)cache[0];
-                        out[x * 4 + 3] = 255;
-                    } else {
-                        out[x] = (UINT8)cache[0];
-                    }
-                }
-            }
-            ImagingSectionLeave(&cookie);
-        }
-        if (inpalette != palette) {
-            ImagingPaletteCacheDelete(palette);
-        }
-    }
-
-    if (inpalette != palette) {
-        ImagingPaletteDelete(palette);
-    }
-
-    return imOut;
-}
 
 static Imaging
 tobilevel(Imaging imOut, Imaging imIn) {
@@ -1368,7 +1169,7 @@ tobilevel(Imaging imOut, Imaging imIn) {
         /* map each pixel to black or white, using error diffusion */
         ImagingSectionEnter(&cookie);
         for (y = 0; y < imIn->ysize; y++) {
-            int l, l0, l1, l2, d2;
+            int l, l0, l1;
             UINT8 *in = (UINT8 *)imIn->image[y];
             UINT8 *out = imOut->image8[y];
 
@@ -1381,14 +1182,10 @@ tobilevel(Imaging imOut, Imaging imIn) {
 
                 /* propagate errors */
                 l -= (int)out[x];
-                l2 = l;
-                d2 = l + l;
-                l += d2;
-                errors[x] = l + l0;
-                l += d2;
-                l0 = l + l1;
-                l1 = l2;
-                l += d2;
+                errors[x] = 3 * l + l0;
+                l0 = 5 * l + l1;
+                l1 = l;
+                l = 7 * l;
             }
 
             errors[x] = l0;
@@ -1399,7 +1196,7 @@ tobilevel(Imaging imOut, Imaging imIn) {
         /* map each pixel to black or white, using error diffusion */
         ImagingSectionEnter(&cookie);
         for (y = 0; y < imIn->ysize; y++) {
-            int l, l0, l1, l2, d2;
+            int l, l0, l1;
             UINT8 *in = (UINT8 *)imIn->image[y];
             UINT8 *out = imOut->image8[y];
 
@@ -1412,14 +1209,10 @@ tobilevel(Imaging imOut, Imaging imIn) {
 
                 /* propagate errors */
                 l -= (int)out[x];
-                l2 = l;
-                d2 = l + l;
-                l += d2;
-                errors[x] = l + l0;
-                l += d2;
-                l0 = l + l1;
-                l1 = l2;
-                l += d2;
+                errors[x] = 3 * l + l0;
+                l0 = 5 * l + l1;
+                l1 = l;
+                l = 7 * l;
             }
 
             errors[x] = l0;