def deltaE_ciede2000(lab1, lab2, kL=1, kC=1, kH=1):
    L1, a1, b1 = np.rollaxis(lab1, -1)[:3]
    L2, a2, b2 = np.rollaxis(lab2, -1)[:3]

    c1 = np.sqrt(a1 ** 2 + b1 ** 2)
    c2 = np.sqrt(a2 ** 2 + b2 ** 2)
    cbar = 0.5 * (c1 + c2)
    c7 = cbar ** 7
    G = 0.5 * (1 - np.sqrt(c7 / (c7 + 25 ** 7)))

    dL_prime = L2 - L1
    Lbar = 0.5 * (L1 + L2)

    a1_prime = a1 * (1 + G)
    a2_prime = a2 * (1 + G)

    c1_prime = np.sqrt(a1_prime ** 2 + b1 ** 2)
    c2_prime = np.sqrt(a2_prime ** 2 + b2 ** 2)
    cbar_prime = 0.5 * (c1_prime + c2_prime)
    dC_prime = c2_prime - c1_prime

    h1_prime = _arctan2pi(b1, a1_prime)
    h2_prime = _arctan2pi(b2, a2_prime)

    dh_prime = h2_prime - h1_prime

    cc = c1_prime * c2_prime
    mask1 = cc == 0.
    mask2 = np.logical_and(-mask1, dh_prime > np.pi)
    mask3 = np.logical_and(-mask1, dh_prime < -np.pi)
    dh_prime = np.where(mask1, 0., dh_prime)
    dh_prime += np.where(mask2, 2 * np.pi, 0)
    dh_prime -= np.where(mask3, 2 * np.pi, 0)

    dH_prime = 2 * np.sqrt(cc) * np.sin(dh_prime / 2)

    Hbar_prime = h1_prime + h2_prime
    mask0 = np.logical_and(np.abs(h1_prime - h2_prime) > np.pi, cc != 0.)
    mask1 = np.logical_and(mask0, Hbar_prime < 2 * np.pi)
    mask2 = np.logical_and(mask0, Hbar_prime >= 2 * np.pi)

    Hbar_prime += np.where(mask1, 2 * np.pi, 0)
    Hbar_prime -= np.where(mask2, 2 * np.pi, 0)
    Hbar_prime *= np.where(cc == 0., 2, 1)
    Hbar_prime *= 0.5

    T = (1 -
         0.17 * np.cos(Hbar_prime - 30 * DEG) +
         0.24 * np.cos(2 * Hbar_prime) +
         0.32 * np.cos(3 * Hbar_prime + 6 * DEG) -
         0.20 * np.cos(4 * Hbar_prime - 63 * DEG)
         )
    dTheta = 30 * DEG * np.exp(-((Hbar_prime / DEG - 275) / 25) ** 2)
    c7 = cbar_prime ** 7
    Rc = 2 * np.sqrt(c7 / (c7 + 25 ** 7))

    term = (Lbar - 50) ** 2
    SL = 1 + 0.015 * term / np.sqrt(20 + term)
    SC = 1 + 0.045 * cbar_prime
    SH = 1 + 0.015 * cbar_prime * T

    RT = -np.sin(2 * dTheta) * Rc

    l_term = dL_prime / (kL * SL)
    c_term = dC_prime / (kC * SC)
    h_term = dH_prime / (kH * SH)
    r_term = RT * c_term * h_term

    dE2 = l_term ** 2
    dE2 += c_term ** 2
    dE2 += h_term ** 2
    dE2 += r_term
    return np.sqrt(dE2)
Color difference as given by the CIEDE 2000 standard.
dataset/ETHPy150Open mrterry/yoink/yoink/delta_e.py/deltaE_ciede2000