diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py index 75f4eba2..b4b9b4c6 100644 --- a/skimage/color/colorconv.py +++ b/skimage/color/colorconv.py @@ -1075,13 +1075,20 @@ def lab2lch(lab): lch = _prepare_lab_array(lab) a, b = lch[..., 1], lch[..., 2] - lch[..., 1], lch[..., 2] = np.sqrt(a ** 2 + b ** 2), np.arctan2(b, a) - - H = lch[..., 2] - H += np.where(H < 0, 2*np.pi, 0) # (-pi, pi) -> (0, 2*pi) + lch[..., 1], lch[..., 2] = _cart2polar_2pi(a, b) return lch +def _cart2polar_2pi(x, y): + """convert cartesian coordiantes to polar (uses non-standard theta range!) + + NON-STANDARD RANGE! Maps to (0, 2*pi) rather than usual (-pi, +pi) + """ + r, t = np.hypot(x, y), np.arctan2(y, x) + t += np.where(t < 0., 2 * np.pi, 0) + return r, t + + def lch2lab(lch): """CIE-LCH to CIE-LAB color space conversion. diff --git a/skimage/color/delta_e.py b/skimage/color/delta_e.py index 2c9d2b19..602e8b66 100644 --- a/skimage/color/delta_e.py +++ b/skimage/color/delta_e.py @@ -21,12 +21,7 @@ from __future__ import division import numpy as np - -def _arctan2pi(b, a): - """np.arctan2 mapped to (0, 2 * pi)""" - ans = np.arctan2(b, a) - ans += np.where(ans < 0, 2 * np.pi, 0.) - return ans +from skimage.color.colorconv import lab2lch, _cart2polar_2pi def deltaE_cie76(lab1, lab2): @@ -50,9 +45,9 @@ def deltaE_cie76(lab1, lab2): .. [2] A. R. Robertson, "The CIE 1976 color-difference formulae," Color Res. Appl. 2, 7-11 (1977). """ - l1, a1, b1 = np.rollaxis(lab1, -1)[:3] - l2, a2, b2 = np.rollaxis(lab2, -1)[:3] - return np.sqrt((l2 - l1) ** 2 + (a2 - a1) ** 2 + (b2 - b1) ** 2) + L1, a1, b1 = np.rollaxis(lab1, -1)[:3] + L2, a2, b2 = np.rollaxis(lab2, -1)[:3] + return np.sqrt((L2 - L1) ** 2 + (a2 - a1) ** 2 + (b2 - b1) ** 2) def deltaE_ciede94(lab1, lab2, kH=1, kC=1, kL=1, k1=0.045, k2=0.015): @@ -104,22 +99,20 @@ def deltaE_ciede94(lab1, lab2, kH=1, kC=1, kL=1, k1=0.045, k2=0.015): .. [1] http://en.wikipedia.org/wiki/Color_difference .. [2] http://www.brucelindbloom.com/index.html?Eqn_DeltaE_CIE94.html """ - l1, a1, b1 = np.rollaxis(lab1, -1)[:3] - l2, a2, b2 = np.rollaxis(lab2, -1)[:3] + L1, C1, h1 = np.rollaxis(lab2lch(lab1), -1)[:3] + L2, C2, h2 = np.rollaxis(lab2lch(lab2), -1)[:3] - dl = l1 - l2 - c1 = np.hypot(a1, b1) - c2 = np.hypot(a2, b2) - dc = c1 - c2 - dh_ab = np.sqrt(deltaE_cie76(lab1, lab2) ** 2 - dl ** 2 - dc ** 2) + dL = L1 - L2 + dC = C1 - C2 + dH_ab = np.sqrt(deltaE_cie76(lab1, lab2) ** 2 - dL ** 2 - dC ** 2) SL = 1 - SC = 1 + k1 * c1 - SH = 1 + k2 * c1 + SC = 1 + k1 * C1 + SH = 1 + k2 * C1 - ans = ((dl / (kL * SL)) ** 2 + - (dc / (kC * SC)) ** 2 + - (dh_ab / (kH * SH)) ** 2 + ans = ((dL / (kL * SL)) ** 2 + + (dC / (kC * SC)) ** 2 + + (dH_ab / (kH * SH)) ** 2 ) return np.sqrt(ans) @@ -166,76 +159,74 @@ 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.hypot(a1, b1) - c2 = np.hypot(a2, b2) - cbar = 0.5 * (c1 + c2) - c7 = cbar ** 7 + # distort `a` based on average chroma + # then convert to lch coordines from distorted `a` + # all subsequence calculations are in the new coordiantes + # (often denoted "prime" in the literature) + Cbar = 0.5 * (np.hypot(a1, b1) + np.hypot(a2, b2)) + c7 = Cbar ** 7 G = 0.5 * (1 - np.sqrt(c7 / (c7 + 25 ** 7))) + scale = 1 + G + C1, h1 = _cart2polar_2pi(a1 * scale, b1) + C2, h2 = _cart2polar_2pi(a2 * scale, b2) + # recall that c, h are polar coordiantes. c==r, h==theta - dL_prime = L2 - L1 + # cide2000 has four terms to delta_e: + # 1) Luminance term + # 2) Hue term + # 3) Chroma term + # 4) hue Rotation term + + # luminance term Lbar = 0.5 * (L1 + L2) + tmp = (Lbar - 50) ** 2 + SL = 1 + 0.015 * tmp / np.sqrt(20 + tmp) + L_term = (L2 - L1) / (kL * SL) - a1_prime = a1 * (1 + G) - a2_prime = a2 * (1 + G) + # chroma term + Cbar = 0.5 * (C1 + C2) # new coordiantes + SC = 1 + 0.045 * Cbar + C_term = (C2 - C1) / (kC * SC) - c1_prime = np.hypot(a1_prime, b1) - c2_prime = np.hypot(a2_prime, b2) - cbar_prime = 0.5 * (c1_prime + c2_prime) - dC_prime = c2_prime - c1_prime + # hue term + h_diff = h2 - h1 + h_sum = h1 + h2 + CC = C1 * C2 - h1_prime = _arctan2pi(b1, a1_prime) - h2_prime = _arctan2pi(b2, a2_prime) + dH = h_diff.copy() + dH[h_diff > np.pi] -= 2 * np.pi + dH[h_diff < -np.pi] += 2 * np.pi + dH[CC == 0.] = 0. # if r == 0, dtheta == 0 + dH_term = 2 * np.sqrt(CC) * np.sin(dH / 2) - 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 + Hbar = h_sum.copy() + mask = np.logical_and(CC != 0., np.abs(h_diff) > np.pi) + Hbar[mask * (h_sum < 2 * np.pi)] += 2 * np.pi + Hbar[mask * (h_sum >= 2 * np.pi)] -= 2 * np.pi + Hbar[CC == 0.] *= 2 + Hbar *= 0.5 T = (1 - - 0.17 * np.cos(Hbar_prime - np.deg2rad(30)) + - 0.24 * np.cos(2 * Hbar_prime) + - 0.32 * np.cos(3 * Hbar_prime + np.deg2rad(6)) - - 0.20 * np.cos(4 * Hbar_prime - np.deg2rad(63)) + 0.17 * np.cos(Hbar - np.deg2rad(30)) + + 0.24 * np.cos(2 * Hbar) + + 0.32 * np.cos(3 * Hbar + np.deg2rad(6)) - + 0.20 * np.cos(4 * Hbar - np.deg2rad(63)) ) - dTheta = (np.deg2rad(30) * - np.exp(-((np.rad2deg(Hbar_prime) - 275) / 25) ** 2) - ) - c7 = cbar_prime ** 7 + SH = 1 + 0.015 * Cbar * T + + H_term = dH_term / (kH * SH) + + # hue rotation + c7 = Cbar ** 7 Rc = 2 * np.sqrt(c7 / (c7 + 25 ** 7)) + dtheta = np.deg2rad(30) * np.exp(-((np.rad2deg(Hbar) - 275) / 25) ** 2) + R_term = -np.sin(2 * dtheta) * Rc * C_term * H_term - 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 + # put it all together + dE2 = L_term ** 2 + dE2 += C_term ** 2 + dE2 += H_term ** 2 + dE2 += R_term return np.sqrt(dE2) @@ -277,27 +268,22 @@ def deltaE_cmc(lab1, lab2, kL=1, kC=1): JPC79 colour-difference formula," J. Soc. Dyers Colour. 100, 128-132 (1984). """ - l1, a1, b1 = np.rollaxis(lab1, -1)[:3] - l2, a2, b2 = np.rollaxis(lab2, -1)[:3] + L1, C1, h1 = np.rollaxis(lab2lch(lab1), -1)[:3] + L2, C2, h2 = np.rollaxis(lab2lch(lab2), -1)[:3] - c1 = np.hypot(a1, b1) - c2 = np.hypot(a2, b2) - dC = c1 - c2 - dl = l1 - l2 - dH = np.sqrt(deltaE_cie76(lab1, lab2) ** 2 - dl ** 2 - dC ** 2) + dC = C1 - C2 + dL = L1 - L2 + dH = np.sqrt(deltaE_cie76(lab1, lab2) ** 2 - dL ** 2 - dC ** 2) - dL = l1 - l2 - - h1 = _arctan2pi(b1, a1) T = np.where(np.logical_and(np.rad2deg(h1) >= 164, np.rad2deg(h1) <= 345), 0.56 + 0.2 * np.abs(np.cos(h1 + np.deg2rad(168))), 0.36 + 0.4 * np.abs(np.cos(h1 + np.deg2rad(35))) ) - c1_4 = c1 ** 4 + c1_4 = C1 ** 4 F = np.sqrt(c1_4 / (c1_4 + 1900)) - SL = np.where(l1 < 16, 0.511, 0.040975 * l1 / (1. + 0.01765 * l1)) - SC = 0.638 + 0.0638 * c1 / (1. + 0.0131 * c1) + SL = np.where(L1 < 16, 0.511, 0.040975 * L1 / (1. + 0.01765 * L1)) + SC = 0.638 + 0.0638 * C1 / (1. + 0.0131 * C1) SH = SC * (F * T + 1 - F) dE2 = (dL / (kL * SL)) ** 2