diff --git a/skimage/color/__init__.py b/skimage/color/__init__.py index 19620cb1..2d072808 100644 --- a/skimage/color/__init__.py +++ b/skimage/color/__init__.py @@ -44,6 +44,12 @@ from .colorconv import (convert_colorspace, from .colorlabel import color_dict, label2rgb +from .delta_e import (deltaE_cie76, + deltaE_ciede94, + deltaE_ciede2000, + deltaE_cmc, + ) + __all__ = ['convert_colorspace', 'guess_spatial_dimensions', @@ -89,4 +95,6 @@ __all__ = ['convert_colorspace', 'is_rgb', 'is_gray', 'color_dict', - 'label2rgb'] + 'label2rgb', + 'deltaE_ciede2000', + ] diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py index d2b20316..36fbeb1a 100644 --- a/skimage/color/colorconv.py +++ b/skimage/color/colorconv.py @@ -26,10 +26,13 @@ Supported color spaces Derived from the RGB CIE color space. Chosen such that ``x == y == z == 1/3`` at the whitepoint, and all color matching functions are greater than zero everywhere. +* LAB CIE : Lightness, a, b +* LCH CIE : Lightness, Chroma, Hue :author: Nicolas Pinto (rgb2hsv) :author: Ralf Gommers (hsv2rgb) :author: Travis Oliphant (XYZ and RGB CIE functions) +:author: Matt Terry (lab2lch) :license: modified BSD @@ -1026,3 +1029,94 @@ def combine_stains(stains, conv_matrix): logrgb2 = np.dot(-np.reshape(stains, (-1, 3)), conv_matrix) rgb2 = np.exp(logrgb2) return rescale_intensity(np.reshape(rgb2 - 2, stains.shape), in_range=(-1, 1)) + + +def lab2lch(lab): + """CIE-LAB to LCH color space conversion. + TODO: something about cylindrical representation + + Parameters + ---------- + lab : array_like + The image in CIE-LAB format, in a 3-D array of shape (.., .., 3). + + Returns + ------- + out : ndarray + The image in LCH format, in a 3-D array of shape (.., .., 3). + + Raises + ------ + ValueError + If `rgb` is not a 3-D array of shape (.., .., 3). + + References + ---------- + + Notes + ----- + The Hue is expressed as an angle between (0, 2*pi) + + Examples + -------- + >>> from skimage import data + >>> from skimage.color import rgb2lab, lab2lch + >>> lena = data.lena() + >>> lena_lab = rgb2lab(lena) + >>> lena_lch = lab2lch(lena_lab) + """ + shape = lab.shape + if len(shape) != 3 or shape[2] != 3: + raise ValueError("Input image expected to be LAB") + + lab = _prepare_colorarray(lab) + lch = np.empty_like(lab) + + a, b = lab[:, :, 1], lab[:, :, 2] + lch[:, :, 0] = lab[:, :, 0] + lch[:, :, 1] = np.sqrt(a**2 + b**2) # C + + H = lch[:, :, 2] = np.arctan2(b, a) + H[H < 0] += 2*np.pi # (-pi, pi) -> (0, 2*pi) + return lch + + +def lch2lab(lch): + """CIE-LCH to CIE-LAB color space conversion. + TODO: something about cylindrical representation + + Parameters + ---------- + lab : array_like + The image in CIE-LCH format, in a 3-D array of shape (.., .., 3). + + Returns + ------- + out : ndarray + The image in LAB format, in a 3-D array of shape (.., .., 3). + + Raises + ------ + ValueError + If `rgb` is not a 3-D array of shape (.., .., 3). + + References + ---------- + + Examples + -------- + >>> from skimage import data + >>> from skimage.color import rgb2lab, lch2lab + >>> lena = data.lena() + >>> lena_lab = rgb2lab(lena) + >>> lena_lch = lab2lch(lena_lab) + >>> lena_lab2 = lch2lab(lena_lch) + """ + lch = _prepare_colorarray(lch) + lab = np.empty_like(lch) + + c, h = lch[:, :, 1], lch[:, :, 2] + lab[:, :, 0] = lch[:, :, 0] + lab[:, :, 1] = c*np.cos(h) + lab[:, :, 2] = c*np.sin(h) + return lab diff --git a/skimage/color/delta_e.py b/skimage/color/delta_e.py new file mode 100644 index 00000000..91072aff --- /dev/null +++ b/skimage/color/delta_e.py @@ -0,0 +1,215 @@ +""" +Functions for calculating the "distance" between colors. + +Implicit in these definitions of "distance" is the notion of "Just Noticible +Distance" (JND). This represents the distance between colors where a human can +percieve different colors. Humans are more sensitive to certain colors than +others, which different deltaE metrics correct for this with varying degrees of +sophistication. + +:author: Matt Terry + +:license: modified BSD + +Reference +--------- + +""" +from __future__ import division + +import numpy as np +from .colorconv import lab2lch + +DEG = np.pi/180 + + +def _unpack_last(x): + x = np.asarray(x) + shape = x.shape + return [x[..., i] for i in range(shape[-1])] + + +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.) + assert ans.max() <= 2*np.pi + assert ans.min() >= 0. + return ans + + +def deltaE_cie76(lab1, lab2): + """ + "just noticible difference" ~ 2.3 + """ + l1, a1, b1 = _unpack_last(lab1) + l2, a2, b2 = _unpack_last(lab2) + return np.sqrt((l2-l1)**2 + (a2-a1)**2 + (b2-b1)**2) + + +def deltaE_ciede94(lab1, lab2, kC=1, kH=1, kL=1, k1=0.045, k2=0.015): + """ + kC, kH are weighting factors, usually unity (default) + + kL, k1, k2 depend on the application. Sample values are: + kL, k1, k2 = 1, 0.045, 0.015 (graphic arts, default) + kL, k1, k2 = 2, 0.048, 0.014 (textiles) + + Note: deltaE_ciede94 the defines the scales for the lightness, hue, and + chroma in terms of the first color. Consequently + deltaE_ciede94(lab1, lab2) != deltaE_ciede94(lab2, lab1) + """ + l1, a1, b1 = _unpack_last(lab1) + l2, a2, b2 = _unpack_last(lab2) + + dl = l1 - l2 + c1 = np.sqrt(a1**2 + b1**2) + c2 = np.sqrt(a2**2 + b2**2) + da = a1 - a2 + db = b1 - b2 + dc = c1 - c2 + dh_ab = np.sqrt(da**2 + db**2 + dc**2) + + SL = 1 + SC = 1 + k1*c1 + SH = 1 + k2*c1 + + ans = (dl/(kL*SL))**2 + (dc/(kC*SC))**2 + (dh_ab/(kH*SH))**2 + return np.sqrt(ans) + + +def deltaE_ciede2000(lab1, lab2, kL=1, kC=1, kH=1): + """ + Parameters + ---------- + lab1 : array_like + pass + lab2 : array_like + pass + kL : float (range), optional + pass + kC : float (range), optional + pass + kH : float (range), optional + pass + + Returns + ------- + deltaE : array_like + The distance between `lab1` and `lab2` + + Notes + ----- + CIEDE 2000 assumes parametric weighting factors for the luminance, chroma, + and hue (kL, kC, kH respectively). + kL = 1 # graphic arts + kL = 2 # textiles + + References + ---------- + http://www.ece.rochester.edu/~gsharma/ciede2000/ciede2000noteCRNA.pdf + """ + L1, a1, b1 = _unpack_last(lab1) + L2, a2, b2 = _unpack_last(lab2) + + 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 = (-mask1) * (dh_prime > np.pi) + mask3 = (-mask1) * (dh_prime < -np.pi) + dh_prime[mask1] = 0. + dh_prime[mask2] += 2*np.pi + dh_prime[mask3] -= 2*np.pi + + 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[mask1] += 2*np.pi + Hbar_prime[mask2] -= 2*np.pi + Hbar_prime[cc == 0.] *= 2 + Hbar_prime *= 0.5 + + deg = np.pi/180. + 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) + + +def deltaE_cmc(lab1, lab2): + """ + indistinguishable if < 1 + usual value for "different" is > 2 + + Note: deltaE_cmc the defines the scales for the lightness, hue, and chroma + in terms of the first color. Consequently + deltaE_cmc(lab1, lab2) != deltaE_cmc(lab2, lab1) + """ + l1, c1, h1 = _unpack_last(lab2lch(lab1)) + l2, c2, h2 = _unpack_last(lab2lch(lab2)) + + sl = np.where(l1 < 16, 0.511, 0.040975*l1 / (1 + 0.01765*l1)) + sc = 0.638 + 0.0638*c1 / (1 + 0.0131*c1) + + c1_4 = c1**4 + f = np.sqrt(c1_4 / (c1_4 + 1900)) + t = np.where(np.logical_and(h1 >= 2.862, h1 <= 6.021), + 0.56 * 0.2 * np.abs(np.cos(h1 + 2.93)), + 0.36 + 0.4 * np.abs(np.cos(h1 + 0.611)) + ) + sh = sc * (f*t + 1-f) + + l, c = 1, 1 + ans = ((l2 - l1)/(l * sl))**2 + ans += ((c2 - c1)/(c * sc))**2 + deg = np.pi/180. + ans += ((h2 - h1)*deg/sh)**2 # metric defines h in terms of degrees + + return np.sqrt(ans) diff --git a/skimage/color/tests/ciede2000_test_data.txt b/skimage/color/tests/ciede2000_test_data.txt new file mode 100644 index 00000000..28503ce5 --- /dev/null +++ b/skimage/color/tests/ciede2000_test_data.txt @@ -0,0 +1,37 @@ +# input, intermediate, and output values for CIEDE2000 dE function +# data taken from "The CIEDE2000 Color-Difference Formula: Implementation Notes, ..." http://www.ece.rochester.edu/~gsharma/ciede2000/ciede2000noteCRNA.pdf +# pair 1 L1 a1 b1 ap1 cp1 hp1 hbar1 G T SL SC SH RT dE 2 L2 a2 b2 ap2 cp2 hp2 +1 1 50.0000 2.6772 -79.7751 2.6774 79.8200 271.9222 270.9611 0.0001 0.6907 1.0000 4.6578 1.8421 -1.7042 2.0425 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000 +2 1 50.0000 3.1571 -77.2803 3.1573 77.3448 272.3395 271.1698 0.0001 0.6843 1.0000 4.6021 1.8216 -1.7070 2.8615 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000 +3 1 50.0000 2.8361 -74.0200 2.8363 74.0743 272.1944 271.0972 0.0001 0.6865 1.0000 4.5285 1.8074 -1.7060 3.4412 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000 +4 1 50.0000 -1.3802 -84.2814 -1.3803 84.2927 269.0618 269.5309 0.0001 0.7357 1.0000 4.7584 1.9217 -1.6809 1.0000 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000 +5 1 50.0000 -1.1848 -84.8006 -1.1849 84.8089 269.1995 269.5997 0.0001 0.7335 1.0000 4.7700 1.9218 -1.6822 1.0000 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000 +6 1 50.0000 -0.9009 -85.5211 -0.9009 85.5258 269.3964 269.6982 0.0001 0.7303 1.0000 4.7862 1.9217 -1.6840 1.0000 2 50.0000 0.0000 -82.7485 0.0000 82.7485 270.0000 +7 1 50.0000 0.0000 0.0000 0.0000 0.0000 0.0000 126.8697 0.5000 1.2200 1.0000 1.0562 1.0229 0.0000 2.3669 2 50.0000 -1.0000 2.0000 -1.5000 2.5000 126.8697 +8 1 50.0000 -1.0000 2.0000 -1.5000 2.5000 126.8697 126.8697 0.5000 1.2200 1.0000 1.0562 1.0229 0.0000 2.3669 2 50.0000 0.0000 0.0000 0.0000 0.0000 0.0000 +9 1 50.0000 2.4900 -0.0010 3.7346 3.7346 359.9847 269.9854 0.4998 0.7212 1.0000 1.1681 1.0404 -0.0022 7.1792 2 50.0000 -2.4900 0.0009 -3.7346 3.7346 179.9862 +10 1 50.0000 2.4900 -0.0010 3.7346 3.7346 359.9847 269.9847 0.4998 0.7212 1.0000 1.1681 1.0404 -0.0022 7.1792 2 50.0000 -2.4900 0.0010 -3.7346 3.7346 179.9847 +11 1 50.0000 2.4900 -0.0010 3.7346 3.7346 359.9847 89.9839 0.4998 0.6175 1.0000 1.1681 1.0346 0.0000 7.2195 2 50.0000 -2.4900 0.0011 -3.7346 3.7346 179.9831 +12 1 50.0000 2.4900 -0.0010 3.7346 3.7346 359.9847 89.9831 0.4998 0.6175 1.0000 1.1681 1.0346 0.0000 7.2195 2 50.0000 -2.4900 0.0012 -3.7346 3.7346 179.9816 +13 1 50.0000 -0.0010 2.4900 -0.0015 2.4900 90.0345 180.0328 0.4998 0.9779 1.0000 1.1121 1.0365 0.0000 4.8045 2 50.0000 0.0009 -2.4900 0.0013 2.4900 270.0311 +14 1 50.0000 -0.0010 2.4900 -0.0015 2.4900 90.0345 180.0345 0.4998 0.9779 1.0000 1.1121 1.0365 0.0000 4.8045 2 50.0000 0.0010 -2.4900 0.0015 2.4900 270.0345 +15 1 50.0000 -0.0010 2.4900 -0.0015 2.4900 90.0345 0.0362 0.4998 1.3197 1.0000 1.1121 1.0493 0.0000 4.7461 2 50.0000 0.0011 -2.4900 0.0016 2.4900 270.0380 +16 1 50.0000 2.5000 0.0000 3.7496 3.7496 0.0000 315.0000 0.4998 0.8454 1.0000 1.1406 1.0396 -0.0001 4.3065 2 50.0000 0.0000 -2.5000 0.0000 2.5000 270.0000 +17 1 50.0000 2.5000 0.0000 3.4569 3.4569 0.0000 346.2470 0.3827 1.4453 1.1608 1.9547 1.4599 -0.0003 27.1492 2 73.0000 25.0000 -18.0000 34.5687 38.9743 332.4939 +18 1 50.0000 2.5000 0.0000 3.4954 3.4954 0.0000 51.7766 0.3981 0.6447 1.0640 1.7498 1.1612 0.0000 22.8977 2 61.0000 -5.0000 29.0000 -6.9907 29.8307 103.5532 +19 1 50.0000 2.5000 0.0000 3.5514 3.5514 0.0000 272.2362 0.4206 0.6521 1.0251 1.9455 1.2055 -0.8219 31.9030 2 56.0000 -27.0000 -3.0000 -38.3556 38.4728 184.4723 +20 1 50.0000 2.5000 0.0000 3.5244 3.5244 0.0000 11.9548 0.4098 1.1031 1.0400 1.9120 1.3353 0.0000 19.4535 2 58.0000 24.0000 15.0000 33.8342 37.0102 23.9095 +21 1 50.0000 2.5000 0.0000 3.7494 3.7494 0.0000 3.5056 0.4997 1.2616 1.0000 1.1923 1.0808 0.0000 1.0000 2 50.0000 3.1736 0.5854 4.7596 4.7954 7.0113 +22 1 50.0000 2.5000 0.0000 3.7493 3.7493 0.0000 0.0000 0.4997 1.3202 1.0000 1.1956 1.0861 0.0000 1.0000 2 50.0000 3.2972 0.0000 4.9450 4.9450 0.0000 +23 1 50.0000 2.5000 0.0000 3.7497 3.7497 0.0000 5.8190 0.4999 1.2197 1.0000 1.1486 1.0604 0.0000 1.0000 2 50.0000 1.8634 0.5757 2.7949 2.8536 11.6380 +24 1 50.0000 2.5000 0.0000 3.7493 3.7493 0.0000 1.9603 0.4997 1.2883 1.0000 1.1946 1.0836 0.0000 1.0000 2 50.0000 3.2592 0.3350 4.8879 4.8994 3.9206 +25 1 60.2574 -34.0099 36.2677 -34.0678 49.7590 133.2085 132.0835 0.0017 1.3010 1.1427 3.2946 1.9951 0.0000 1.2644 2 60.4626 -34.1751 39.4387 -34.2333 52.2238 130.9584 +26 1 63.0109 -31.0961 -5.8663 -32.6194 33.1427 190.1951 188.8221 0.0490 0.9402 1.1831 2.4549 1.4560 0.0000 1.2630 2 62.8187 -29.7946 -4.0864 -31.2542 31.5202 187.4490 +27 1 61.2901 3.7196 -5.3901 5.5668 7.7487 315.9240 310.0313 0.4966 0.6952 1.1586 1.3092 1.0717 -0.0032 1.8731 2 61.4292 2.2480 -4.9620 3.3644 5.9950 304.1385 +28 1 35.0831 -44.1164 3.7933 -44.3939 44.5557 175.1161 176.4290 0.0063 1.0168 1.2148 2.9105 1.6476 0.0000 1.8645 2 35.0232 -40.0716 1.5901 -40.3237 40.3550 177.7418 +29 1 22.7233 20.0904 -46.6940 20.1424 50.8532 293.3339 291.3809 0.0026 0.3636 1.4014 3.1597 1.2617 -1.2537 2.0373 2 23.0331 14.9730 -42.5619 15.0118 45.1317 289.4279 +30 1 36.4612 47.8580 18.3852 47.9197 51.3256 20.9901 21.8781 0.0013 0.9239 1.1943 3.3888 1.7357 0.0000 1.4146 2 36.2715 50.5065 21.2231 50.5716 54.8444 22.7660 +31 1 90.8027 -2.0831 1.4410 -3.1245 3.4408 155.2410 167.1011 0.4999 1.1546 1.6110 1.1329 1.0511 0.0000 1.4441 2 91.1528 -1.6435 0.0447 -2.4651 2.4655 178.9612 +32 1 90.9257 -0.5406 -0.9208 -0.8109 1.2270 228.6315 218.4363 0.5000 1.3916 1.5930 1.0620 1.0288 0.0000 1.5381 2 88.6381 -0.8985 -0.7239 -1.3477 1.5298 208.2412 +33 1 6.7747 -0.2908 -2.4247 -0.4362 2.4636 259.8025 263.0049 0.4999 0.9556 1.6517 1.1057 1.0337 -0.0004 0.6377 2 5.8714 -0.0985 -2.2286 -0.1477 2.2335 266.2073 +34 1 2.0776 0.0795 -1.1350 0.1192 1.1412 275.9978 268.0910 0.5000 0.7826 1.7246 1.0383 1.0100 0.0000 0.9082 2 0.9033 -0.0636 -0.5514 -0.0954 0.5596 260.18421 diff --git a/skimage/color/tests/test_ciede2000.py b/skimage/color/tests/test_ciede2000.py new file mode 100644 index 00000000..5c00fb24 --- /dev/null +++ b/skimage/color/tests/test_ciede2000.py @@ -0,0 +1,61 @@ +"""Test for correctness of color distance functions + +Authors +------- +Matt Terry + +:license: modified BSD +""" +from os.path import abspath, dirname, join as pjoin + +import numpy as np +from numpy.testing import assert_array_almost_equal + +from skimage.color import deltaE_ciede2000 + + +def test_ciede2000_dE(): + dtype = [('pair', int), + ('1', int), + ('L1', float), + ('a1', float), + ('b1', float), + ('a1_prime', float), + ('C1_prime', float), + ('h1_prime', float), + ('hbar_prime', float), + ('G', float), + ('T', float), + ('SL', float), + ('SC', float), + ('SH', float), + ('RT', float), + ('dE', float), + ('2', int), + ('L2', float), + ('a2', float), + ('b2', float), + ('a2_prime', float), + ('C2_prime', float), + ('h2_prime', float), + ] + + # note: ciede_test_data.txt contains several intermediate quantities + path = pjoin(dirname(abspath(__file__)), 'ciede_test_data.txt') + data = np.loadtxt(path, dtype=dtype) + + N = len(data) + + lab1 = np.zeros((N, 3)) + lab1[:, 0] = data['L1'] + lab1[:, 1] = data['a1'] + lab1[:, 2] = data['b1'] + + lab2 = np.zeros((N, 3)) + lab2[:, 0] = data['L2'] + lab2[:, 1] = data['a2'] + lab2[:, 2] = data['b2'] + + dE2 = deltaE_ciede2000(lab1, lab2) + + assert_array_almost_equal(dE2, data['dE'])