1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
function dE2000 = ComputeDE2000_Lab(Lab1,Lab2, kLCH)
% Computes the deltaE values (CIEDE2000) for color pairs given in CIELAB coordinates (L*a*b*).
%
% dE2000 = ComputeDE2000_Lab(Lab1,Lab2[,kLCH]);
%
% The inputs 'Lab1' and 'Lab2' are N*3 matrices each with one color per row [L*,a*,b*].
% The output is a column vector with N elements, containing the according
% dE2000 values.
% 'kLCH' is optional and contains one or three elements for the weighting
% factors kL, kC, kH (either one value for all, or three values for each).
% The weighting factors are set to 1, if 'kLCH' is not specified.
%
% This implementation might not be the most efficient but it follows the formulas
% given in https://en.wikipedia.org/wiki/Color_difference#CIEDE2000 as closely
% as possible for better traceability.
% HISTORY:
% 12/01/2017 (MR = Marc Repnow, EPFL) Initial version.
if nargin==0
testMe();
return
end
if any(size(Lab1)~=size(Lab2)) || size(Lab1,2)~=3
error('L*a*b* input parameters must be of size Nx3 (N>=1).');
end
if nargin>2
if length(kLCH)==1
kLCH = kLCH*[1,1,1];
elseif isempty(kLCH)
kLCH = [1,1,1];
elseif numel(kLCH)~=3
error('Input parameter ''kLCH'' must have either 1 or 3 elements.');
end
else
kLCH = [1,1,1];
end
%--- "Degree"-versions of sin() and cos().
sinD = @(deg) sin(deg*pi()/180);
cosD = @(deg) cos(deg*pi()/180);
%---
Lstar12 = [Lab1(:,1), Lab2(:,1)];
aStar12 = [Lab1(:,2), Lab2(:,2)];
bStar12 = [Lab1(:,3), Lab2(:,3)];
CprimeStar12 = sqrt(aStar12.^2 + bStar12.^2);
dLprime = Lstar12(:,2) - Lstar12(:,1);
Lbar = (Lstar12(:,2) + Lstar12(:,1))/2;
Cbar = (CprimeStar12(:,1)+CprimeStar12(:,2))/2;
aPrime12 = aStar12 .* repmat(1 + 0.5*(1 - sqrt(Cbar.^7./(Cbar.^7 + 25^7))), 1,2);
Cprime12 = sqrt(aPrime12.^2 + bStar12.^2);
dCprime = Cprime12(:,2)-Cprime12(:,1);
CbarPrime = (Cprime12(:,1)+Cprime12(:,2))/2;
%--- Note that atan2(0,0) will return 0.
hPrime12 = mod(atan2(bStar12, aPrime12)*180/pi() + 10*360, 360);
N = length(dLprime);
dhPrime = zeros(N,1);
HbarPrime = zeros(N,1);
for i=1:N
diff_hPrime12 = hPrime12(i,2)-hPrime12(i,1);
sum_hPrime12 = hPrime12(i,2)+hPrime12(i,1);
if any(abs(Cprime12(i,:))<10*eps)
dhPrime(i) = 0;
elseif abs(diff_hPrime12)<=180
dhPrime(i) = diff_hPrime12;
elseif diff_hPrime12 <= 0
dhPrime(i) = diff_hPrime12 + 360;
else
dhPrime(i) = diff_hPrime12 - 360;
end
if any(abs(Cprime12(i,:))<10*eps)
HbarPrime(i) = sum_hPrime12;
elseif abs(diff_hPrime12)<=180
HbarPrime(i) = (sum_hPrime12)/2;
elseif sum_hPrime12 < 360
HbarPrime(i) = (sum_hPrime12 + 360)/2;
else
HbarPrime(i) = (sum_hPrime12 - 360)/2;
end
end
dHprime = 2*sqrt(Cprime12(:,1).*Cprime12(:,2)).*sinD(dhPrime/2);
T = 1 - 0.17*cosD(HbarPrime-30) + 0.24*cosD(2*HbarPrime) ...
+ 0.32*cosD(3*HbarPrime+6) - 0.2*cosD(4*HbarPrime-63);
SL = 1 + 0.015*(Lbar-50).^2 ./ sqrt(20 + (Lbar-50).^2);
SC = 1 + 0.045*CbarPrime;
SH = 1 + 0.015*CbarPrime.*T;
RT = -2*sqrt(CbarPrime.^7./(CbarPrime.^7 + 25^7)) .* sinD(60*exp(-((HbarPrime-275)./25).^2));
kL = kLCH(1);
kC = kLCH(2);
kH = kLCH(3);
dE2000 = sqrt((dLprime./(kL*SL)).^2 + (dCprime./(kC*SC)).^2 + (dHprime./(kH*SH)).^2 ...
+ RT.*dCprime./(kC*SC).*dHprime./(kH*SH));
end %%%%
%--------------------------------------------------------------------------------------
function testMe()
% Test data are from
% Sharma et al., 2005, "The CIEDE2000 Color-Difference Formula: Implementation Notes,
% Supplementary Test Data, and Mathematical Observations".
% The values in this publication are given with 4 decimal places. For this test function here,
% the MATLAB script provided by Sharma was used to get more precise values (8 decimal places).
% L*a*b* 1, L*a*b* 2, dE2000
testset = ...
[ 50.0000, 2.6772, -79.7751, 50.0000, 0.0000, -82.7485, 2.04245968;
50.0000, 3.1571, -77.2803, 50.0000, 0.0000, -82.7485, 2.86151017;
50.0000, 2.8361, -74.0200, 50.0000, 0.0000, -82.7485, 3.44119060;
50.0000, -1.3802, -84.2814, 50.0000, 0.0000, -82.7485, 0.99999886;
50.0000, -1.1848, -84.8006, 50.0000, 0.0000, -82.7485, 1.00000470;
50.0000, -0.9009, -85.5211, 50.0000, 0.0000, -82.7485, 1.00001297;
50.0000, 0.0000, 0.0000, 50.0000, -1.0000, 2.0000, 2.36685882;
50.0000, -1.0000, 2.0000, 50.0000, 0.0000, 0.0000, 2.36685882;
50.0000, 2.4900, -0.0010, 50.0000, -2.4900, 0.0009, 7.17917201;
50.0000, 2.4900, -0.0010, 50.0000, -2.4900, 0.0010, 7.17916264;
50.0000, 2.4900, -0.0010, 50.0000, -2.4900, 0.0011, 7.21947215;
50.0000, 2.4900, -0.0010, 50.0000, -2.4900, 0.0012, 7.21947421;
50.0000, -0.0010, 2.4900, 50.0000, 0.0009, -2.4900, 4.80452169;
50.0000, -0.0010, 2.4900, 50.0000, 0.0010, -2.4900, 4.80452451;
50.0000, -0.0010, 2.4900, 50.0000, 0.0011, -2.4900, 4.74607111;
50.0000, 2.5000, 0.0000, 50.0000, 0.0000, -2.5000, 4.30648210;
50.0000, 2.5000, 0.0000, 73.0000, 25.0000, -18.0000, 27.14923130;
50.0000, 2.5000, 0.0000, 61.0000, -5.0000, 29.0000, 22.89769247;
50.0000, 2.5000, 0.0000, 56.0000, -27.0000, -3.0000, 31.90300465;
50.0000, 2.5000, 0.0000, 58.0000, 24.0000, 15.0000, 19.45352143;
50.0000, 2.5000, 0.0000, 50.0000, 3.1736, 0.5854, 1.00002634;
50.0000, 2.5000, 0.0000, 50.0000, 3.2972, 0.0000, 0.99997287;
50.0000, 2.5000, 0.0000, 50.0000, 1.8634, 0.5757, 1.00004950;
50.0000, 2.5000, 0.0000, 50.0000, 3.2592, 0.3350, 1.00003476;
60.2574, -34.0099, 36.2677, 60.4626, -34.1751, 39.4387, 1.26442001;
63.0109, -31.0961, -5.8663, 62.8187, -29.7946, -4.0864, 1.26295930;
61.2901, 3.7196, -5.3901, 61.4292, 2.2480, -4.9620, 1.87307050;
35.0831, -44.1164, 3.7933, 35.0232, -40.0716, 1.5901, 1.86449523;
22.7233, 20.0904, -46.6940, 23.0331, 14.9730, -42.5619, 2.03725827;
36.4612, 47.8580, 18.3852, 36.2715, 50.5065, 21.2231, 1.41457792;
90.8027, -2.0831, 1.4410, 91.1528, -1.6435, 0.0447, 1.44412908;
90.9257, -0.5406, -0.9208, 88.6381, -0.8985, -0.7239, 1.53811701;
6.7747, -0.2908, -2.4247, 5.8714, -0.0985, -2.2286, 0.63772767;
2.0776, 0.0795, -1.1350, 0.9033, -0.0636, -0.5514, 0.90823284];
dE2000 = ComputeDE2000_Lab(testset(:,1:3), testset(:,4:6));
%--- Our test case deltaE's are given with 8 decimal places,
% so ignore differences that can be explained by the according
% round-off errors.
if any(abs(dE2000 - testset(:,end)) > 0.5e-8)
fprintf('ComputeDE2000_Lab test failed.\n');
[~, i] = max(abs(dE2000 - testset(:,end)));
fprintf('e.g. worst test case (#%d): Lab2dE2000([%.4f,%.4f,%.4f],[%.4f,%.4f,%.4f]==%.8f (expected: %.8f)\n',...
i, testset(i,1:6), dE2000(i), testset(i,end));
else
fprintf('ComputeDE2000_Lab test passed.\n');
end
end
|
