Color Similarity + Conversion (Python)

ciede.py

from math import atan2, cos, exp, fabs, hypot, pi, sin, sqrt
from typing import Tuple

lab_type = Tuple[float, float, float]
hex_type = str  # e.g: #ff00ff
rgb_type = Tuple[int, int, int]


def hex2rgb(v: hex_type) -> rgb_type:
    return int(v[1:3], 16), int(v[3:5], 16), int(v[5:7], 16)


def srgb2linear(c: int) -> float:
    c = c / 255.0
    return c / 12.92 if c <= 0.04045 else pow((c + 0.055) / 1.055, 2.4)


def rgb2lab(r: int, g: int, b: int) -> lab_type:
    r, g, b = srgb2linear(r), srgb2linear(g), srgb2linear(b)

    # Convert linear RGB to XYZ (D65)
    X = r * 0.4124564 + g * 0.3575761 + b * 0.1804375
    Y = r * 0.2126729 + g * 0.7151522 + b * 0.0721750
    Z = r * 0.0193339 + g * 0.1191920 + b * 0.9503041

    # Reference white D65
    Xn, Yn, Zn = 0.95047, 1.00000, 1.08883

    def f(t):
        d = 6 / 29
        return pow(t, 1 / 3) if t > pow(d, 3) else t / (3 * d * d) + 4 / 29

    fx, fy, fz = f(X / Xn), f(Y / Yn), f(Z / Zn)

    L = 116 * fy - 16
    a = 500 * (fx - fy)
    b = 200 * (fy - fz)

    return L, a, b


def hex2lab(v: hex_type) -> lab_type:
    return rgb2lab(*hex2rgb(v))


def rad2deg(rad):
    return (360 * rad) / (2 * pi)


def deg2rad(deg):
    return (2 * pi * deg) / 360


def ciede2000(a: lab_type, b: lab_type) -> float:
    """
    The classic CIE ΔE2000 implementation, which operates on two L*a*b* colors,
    and returns their difference.

    credits: https://github.com/michel-leonard

    L1 = 33.0           a1 = 56.551         b1 = 107.34
    L2 = 33.0           a2 = 56.551         b2 = 107.41
    CIE ΔE2000 = ΔE00 = 0.0182199484

    L1 = 57.751         a1 = -64.8572       b1 = 104.6669
    L2 = 57.751         a2 = -65.0          b2 = 104.6669
    CIE ΔE2000 = ΔE00 = 0.04105095889

    L1 = 51.2347        a1 = -53.74         b1 = -85.6
    L2 = 51.2347        a2 = -54.0          b2 = -85.6
    CIE ΔE2000 = ΔE00 = 0.07251289185

    L1 = 27.5           a1 = -38.963        b1 = 102.02
    L2 = 27.5           a2 = -39.0          b2 = 101.6285
    CIE ΔE2000 = ΔE00 = 0.09222135798

    L1 = 58.207         a1 = -105.0         b1 = -79.71
    L2 = 58.207         a2 = -106.2         b2 = -79.71
    CIE ΔE2000 = ΔE00 = 0.23933030702

    L1 = 18.6           a1 = 5.4            b1 = 81.0
    L2 = 18.6           a2 = 5.4            b2 = 90.0
    CIE ΔE2000 = ΔE00 = 1.87770163489

    L1 = 32.913         a1 = 3.0961         b1 = 115.117
    L2 = 32.913         a2 = 8.4101         b2 = 115.117
    CIE ΔE2000 = ΔE00 = 2.60306648026

    L1 = 49.0           a1 = -12.0          b1 = 124.282
    L2 = 53.0           a2 = -12.0          b2 = 116.703
    CIE ΔE2000 = ΔE00 = 4.16877634871

    L1 = 21.032         a1 = -30.4          b1 = -68.783
    L2 = 26.0           a2 = -14.9          b2 = -54.29
    CIE ΔE2000 = ΔE00 = 8.23700340366

    L1 = 76.5           a1 = -83.7311       b1 = 10.0
    L2 = 97.08          a2 = -59.472        b2 = 28.1
    CIE ΔE2000 = ΔE00 = 17.25538447747
    """
    l_1, a_1, b_1 = a
    l_2, a_2, b_2 = b
    # Working in Python with the CIEDE2000 color-difference formula.
    # k_l, k_c, k_h are parametric factors to be adjusted according to
    # different viewing parameters such as textures, backgrounds...
    k_l = k_c = k_h = 1.0
    n = (hypot(a_1, b_1) + hypot(a_2, b_2)) * 0.5
    n = n * n * n * n * n * n * n
    # A factor involving chroma raised to the power of 7 designed to make
    # the influence of chroma on the total color difference more accurate.
    n = 1.0 + 0.5 * (1.0 - sqrt(n / (n + 6103515625.0)))
    # hypot calculates the Euclidean distance while avoiding overflow/underflow
    c_1 = hypot(a_1 * n, b_1)
    c_2 = hypot(a_2 * n, b_2)
    # atan2 is preferred over atan because it accurately computes the angle of
    # a point (x, y) in all quadrants, handling the signs of both coordinates.
    h_1 = atan2(b_1, a_1 * n)
    h_2 = atan2(b_2, a_2 * n)
    h_1 += 2.0 * pi * (h_1 < 0.0)
    h_2 += 2.0 * pi * (h_2 < 0.0)
    n = abs(h_2 - h_1)
    # Cross-implementation consistent rounding.
    if pi - 1e-14 < n and n < pi + 1e-14:
        n = pi
    # When the hue angles lie in different quadrants, the straightforward
    # average can produce a mean that incorrectly suggests a hue angle in
    # the wrong quadrant, the next lines handle this issue.
    h_m = (h_1 + h_2) * 0.5
    h_d = (h_2 - h_1) * 0.5
    if pi < n:
        if 0.0 < h_d:
            h_d -= pi
        else:
            h_d += pi
        h_m += pi
    p = 36.0 * h_m - 55.0 * pi
    n = (c_1 + c_2) * 0.5
    n = n * n * n * n * n * n * n
    # The hue rotation correction term is designed to account for the
    # non-linear behavior of hue differences in the blue region.
    r_t = (
        -2.0
        * sqrt(n / (n + 6103515625.0))
        * sin(pi / 3.0 * exp(p * p / (-25.0 * pi * pi)))
    )
    n = (l_1 + l_2) * 0.5
    n = (n - 50.0) * (n - 50.0)
    # Lightness.
    lt = (l_2 - l_1) / (k_l * (1.0 + 0.015 * n / sqrt(20.0 + n)))
    # These coefficients adjust the impact of different harmonic
    # components on the hue difference calculation.
    t = (
        1.0
        + 0.24 * sin(2.0 * h_m + pi * 0.5)
        + 0.32 * sin(3.0 * h_m + 8.0 * pi / 15.0)
        - 0.17 * sin(h_m + pi / 3.0)
        - 0.20 * sin(4.0 * h_m + 3.0 * pi / 20.0)
    )
    n = c_1 + c_2
    # Hue.
    h = 2.0 * sqrt(c_1 * c_2) * sin(h_d) / (k_h * (1.0 + 0.0075 * n * t))
    # Chroma.
    c = (c_2 - c_1) / (k_c * (1.0 + 0.0225 * n))
    # Returning the square root ensures that the result represents
    # the "true" geometric distance in the color space.
    return sqrt(lt * lt + h * h + c * c + c * h * r_t)

  
def ciede2000_weighted(a: lab_type, b: lab_type, Kl=1, Kc=1, Kh=1) -> float:
    """
    Computes color difference as developed by the International Commission on
    Illumination (CIE) in 2000. The implementation is based on the formula
    from Bruce Lindbloom. Resulting values range from 0 (no difference)
    to 100 (maximum difference), and are a metric for how the human eye
    percieves color difference.
    The optional parameters Kl, Kc, and Kh may be used to adjust
    weightings of lightness, chroma, and hue.

    credits: https://github.com/gka/chroma.js
    """
    L1, a1, b1 = a
    L2, a2, b2 = b

    avgL = (L1 + L2) / 2
    C1 = sqrt(a1 * a1 + b1 * b1)
    C2 = sqrt(a2 * a2 + b2 * b2)
    avgC = (C1 + C2) / 2

    G = 0.5 * (1 - sqrt(pow(avgC, 7) / (pow(avgC, 7) + pow(25, 7))))

    a1p = a1 * (1 + G)
    a2p = a2 * (1 + G)

    C1p = sqrt(a1p * a1p + b1 * b1)
    C2p = sqrt(a2p * a2p + b2 * b2)
    avgCp = (C1p + C2p) / 2

    h1p = rad2deg(atan2(b1, a1p))
    if h1p < 0:
        h1p += 360

    h2p = rad2deg(atan2(b2, a2p))
    if h2p < 0:
        h2p += 360

    avgHp = (h1p + h2p + 360) / 2 if fabs(h1p - h2p) > 180 else (h1p + h2p) / 2

    T = (
        1
        - 0.17 * cos(deg2rad(avgHp - 30))
        + 0.24 * cos(deg2rad(2 * avgHp))
        + 0.32 * cos(deg2rad(3 * avgHp + 6))
        - 0.20 * cos(deg2rad(4 * avgHp - 63))
    )

    deltaHp = h2p - h1p
    if fabs(deltaHp) > 180:
        deltaHp += 360 if h2p <= h1p else -360

    deltaHp = 2 * sqrt(C1p * C2p) * sin(deg2rad(deltaHp / 2))

    deltaL = L2 - L1
    deltaCp = C2p - C1p

    Sl = 1 + (0.015 * pow(avgL - 50, 2)) / sqrt(20 + pow(avgL - 50, 2))
    Sc = 1 + 0.045 * avgCp
    Sh = 1 + 0.015 * avgCp * T

    deltaTheta = 30 * exp(-pow((avgHp - 275) / 25, 2))
    Rc = 2 * sqrt(pow(avgCp, 7) / (pow(avgCp, 7) + pow(25, 7)))
    Rt = -Rc * sin(deg2rad(2 * deltaTheta))

    result = sqrt(
        pow(deltaL / (Kl * Sl), 2)
        + pow(deltaCp / (Kc * Sc), 2)
        + pow(deltaHp / (Kh * Sh), 2)
        + Rt * (deltaCp / (Kc * Sc)) * (deltaHp / (Kh * Sh))
    )

    return max(0, min(100, result))

test_ciede.py

import pytest

from ciede import ciede2000, ciede2000_weighted, hex2lab


# Test cases derived from
# http://www2.ece.rochester.edu/~gsharma/ciede2000/dataNprograms
# /ciede2000testdata.txt
@pytest.mark.parametrize(
    "lab1, lab2, expected",
    [
        ([50.0000, 2.6772, -79.7751], [50.0000, 0.0000, -82.7485], 2.0425),
        ([50.0000, 3.1571, -77.2803], [50.0000, 0.0000, -82.7485], 2.8615),
        ([50.0000, 2.8361, -74.0200], [50.0000, 0.0000, -82.7485], 3.4412),
        ([50.0000, -1.3802, -84.2814], [50.0000, 0.0000, -82.7485], 1.0000),
        ([50.0000, -1.1848, -84.8006], [50.0000, 0.0000, -82.7485], 1.0000),
        ([50.0000, -0.9009, -85.5211], [50.0000, 0.0000, -82.7485], 1.0000),
        ([50.0000, 0.0000, 0.0000], [50.0000, -1.0000, 2.0000], 2.3669),
        ([50.0000, -1.0000, 2.0000], [50.0000, 0.0000, 0.0000], 2.3669),
        ([50.0000, 2.4900, -0.0010], [50.0000, -2.4900, 0.0009], 7.1792),
        ([50.0000, 2.4900, -0.0010], [50.0000, -2.4900, 0.0010], 7.1792),
        ([50.0000, 2.4900, -0.0010], [50.0000, -2.4900, 0.0011], 7.2195),
        ([50.0000, 2.4900, -0.0010], [50.0000, -2.4900, 0.0012], 7.2195),
        ([50.0000, -0.0010, 2.4900], [50.0000, 0.0009, -2.4900], 4.8045),
        ([50.0000, -0.0010, 2.4900], [50.0000, 0.0010, -2.4900], 4.8045),
        ([50.0000, -0.0010, 2.4900], [50.0000, 0.0011, -2.4900], 4.7461),
        ([50.0000, 2.5000, 0.0000], [50.0000, 0.0000, -2.5000], 4.3065),
        ([50.0000, 2.5000, 0.0000], [73.0000, 25.0000, -18.0000], 27.1492),
        ([50.0000, 2.5000, 0.0000], [61.0000, -5.0000, 29.0000], 22.8977),
        ([50.0000, 2.5000, 0.0000], [56.0000, -27.0000, -3.0000], 31.9030),
        ([50.0000, 2.5000, 0.0000], [58.0000, 24.0000, 15.0000], 19.4535),
        ([50.0000, 2.5000, 0.0000], [50.0000, 3.1736, 0.5854], 1.0000),
        ([50.0000, 2.5000, 0.0000], [50.0000, 3.2972, 0.0000], 1.0000),
        ([50.0000, 2.5000, 0.0000], [50.0000, 1.8634, 0.5757], 1.0000),
        ([50.0000, 2.5000, 0.0000], [50.0000, 3.2592, 0.3350], 1.0000),
        ([60.2574, -34.0099, 36.2677], [60.4626, -34.1751, 39.4387], 1.2644),
        ([63.0109, -31.0961, -5.8663], [62.8187, -29.7946, -4.0864], 1.2630),
        ([61.2901, 3.7196, -5.3901], [61.4292, 2.2480, -4.9620], 1.8731),
        ([35.0831, -44.1164, 3.7933], [35.0232, -40.0716, 1.5901], 1.8645),
        ([22.7233, 20.0904, -46.6940], [23.0331, 14.9730, -42.5619], 2.0373),
        ([36.4612, 47.8580, 18.3852], [36.2715, 50.5065, 21.2231], 1.4146),
        ([90.8027, -2.0831, 1.4410], [91.1528, -1.6435, 0.0447], 1.4441),
        ([90.9257, -0.5406, -0.9208], [88.6381, -0.8985, -0.7239], 1.5381),
        ([6.7747, -0.2908, -2.4247], [5.8714, -0.0985, -2.2286], 0.6377),
        ([2.0776, 0.0795, -1.1350], [0.9033, -0.0636, -0.5514], 0.9082),
    ],
)
def test_ciede2000(lab1, lab2, expected):
    result = ciede2000(lab1, lab2)
    assert round(result, 4) == expected

    result = ciede2000_weighted(lab1, lab2)
    assert round(result, 4) == expected


@pytest.mark.parametrize(
    "input_colors, expected",
    (
        ((0x000000, 0x000000), 0),
        ((0xFFFFFF, 0x000000), 100),
        ((0xFF0000, 0x00FF00), 86.6082374535373),
        ((0x00FF00, 0xFF0000), 86.6082374535373),
        ((0x00BEEF, 0x00BEEE), 0.28076037937072745),
        ((0xFFFF00, 0x0000FF), 100),
    ),
)
def test_delta_e(input_colors, expected):
    c1 = hex2lab(f"#{input_colors[0]:06x}")
    c2 = hex2lab(f"#{input_colors[1]:06x}")
    result = ciede2000_weighted(c1, c2)
    assert round(result, 4) == round(expected, 4)