michel-leonard
commented
2 months ago My proposition (avoiding deg/rad conversion) based on your implementation would look like :
def delta_E_CIE2000(Lab1, Lab2):
import numpy as np
# Unpack Lab values
L1, a1, b1 = np.ravel(Lab1)
L2, a2, b2 = np.ravel(Lab2)
# Constants
kL = kC = kH = 1
# Calculate the averages
l_bar_prime = (L1 + L2) / 2
# Calculate C1, C2, and their average
c1 = np.hypot(a1, b1)
c2 = np.hypot(a2, b2)
c_bar = (c1 + c2) / 2
# Calculate G factor
c_bar7 = c_bar ** 7
g = 0.5 * (1 - np.sqrt(c_bar7 / (c_bar7 + 25 ** 7)))
# Calculate a_prime and c_prime values
a1_prime = a1 * (1 + g)
a2_prime = a2 * (1 + g)
c1_prime = np.hypot(a1_prime, b1)
c2_prime = np.hypot(a2_prime, b2)
c_bar_prime = (c1_prime + c2_prime) / 2
# Calculate h_prime values
h1_prime = np.arctan2(b1, a1_prime)
h2_prime = np.arctan2(b2, a2_prime)
# Adjust h_prime values
if h1_prime < 0:
h1_prime += 2 * np.pi
if h2_prime < 0:
h2_prime += 2 * np.pi
# Calculate the average hue
if np.abs(h1_prime - h2_prime) > np.pi:
h_bar_prime = (h1_prime + h2_prime + 2 * np.pi) / 2
else:
h_bar_prime = (h1_prime + h2_prime) / 2
# Calculate T
t = (1
- 0.17 * np.cos(h_bar_prime - np.pi / 6)
+ 0.24 * np.cos(2 * h_bar_prime)
+ 0.32 * np.cos(3 * h_bar_prime + np.pi / 30)
- 0.20 * np.cos(4 * h_bar_prime - 7 * np.pi / 20))
# Calculate delta_h_prime
if np.abs(h2_prime - h1_prime) <= np.pi:
delta_h_prime = h2_prime - h1_prime
elif h2_prime <= h1_prime:
delta_h_prime = h2_prime - h1_prime + 2 * np.pi
else:
delta_h_prime = h2_prime - h1_prime - 2 * np.pi
# Calculate delta values
delta_L_prime = L2 - L1
delta_C_prime = c2_prime - c1_prime
delta_H_prime = 2 * np.sqrt(c1_prime * c2_prime) * np.sin(delta_h_prime / 2)
# Calculate S_L, S_C, S_H
l50_sq = (l_bar_prime - 50) ** 2
sL = 1 + 0.015 * l50_sq / np.sqrt(20 + l50_sq)
sC = 1 + 0.045 * c_bar_prime
sH = 1 + 0.015 * c_bar_prime * t
# Calculate R_C and R_T
c_bar_prime7 = c_bar_prime ** 7
rC = np.sqrt(c_bar_prime7 / (c_bar_prime7 + 25 ** 7))
rT = -2 * rC * np.sin(np.pi / 3 * np.exp(-((36 * h_bar_prime - 55 * np.pi) ** 2) / (25 * np.pi ** 2)))
# Calculate delta_E
delta_E = np.sqrt(
(delta_L_prime / (kL * sL)) ** 2 +
(delta_C_prime / (kC * sC)) ** 2 +
(delta_H_prime / (kH * sH)) ** 2 +
(delta_C_prime / (kC * sC)) * (delta_H_prime / (kH * sH)) * rT
)
return delta_E
KelSolaar
Description
Hello,
I've been using the CIEDE2000 color-difference calculation function (
delta_E_CIE2000) in the library and noticed that the function performs an unnecessary conversion between radians and degrees. Specifically, it multiplies by 180 and then divides by 180 shortly afterward, which seems redundant. This conversion introduces a small epsilon error in certain cases.For example, with the inputs
L1, a1, b1, L2, a2, b2 = 88, -124, 56, 97, 62, -28, the algorithm reaches a condition where it evaluatesnp.fabs(h1_prime - h2_prime)to "180.00000000000003" degrees, which is "π" (3.141592653589793) radians when you omit the conversion from radians to degrees.The specific branch taken right after when the value is greater than π (radians) or 180° is inconsistent due to this minor error, affecting the output of the algorithm.
Thank you for your attention to this detail, and I'm happy to provide further information or assist with testing if needed.
Best regards, Michel Leonard
Environment Information