-
Notifications
You must be signed in to change notification settings - Fork 28
/
Copy pathCIEDE2000.cpp
174 lines (160 loc) · 4.59 KB
/
CIEDE2000.cpp
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
/*
* CIEDE2000.cpp
* Part of http://github.com/gfiumara/CIEDE2000 by Gregory Fiumara.
* See LICENSE for details.
*/
#include <cmath>
#include <CIEDE2000.h>
/*******************************************************************************
* Conversions.
******************************************************************************/
constexpr double
CIEDE2000::deg2Rad(
const double deg)
{
return (deg * (M_PI / 180.0));
}
constexpr double
CIEDE2000::rad2Deg(
const double rad)
{
return ((180.0 / M_PI) * rad);
}
double
CIEDE2000::CIEDE2000(
const LAB &lab1,
const LAB &lab2)
{
/*
* "For these and all other numerical/graphical delta E00 values
* reported in this article, we set the parametric weighting factors
* to unity(i.e., k_L = k_C = k_H = 1.0)." (Page 27).
*/
const double k_L = 1.0, k_C = 1.0, k_H = 1.0;
const double deg360InRad = CIEDE2000::deg2Rad(360.0);
const double deg180InRad = CIEDE2000::deg2Rad(180.0);
const double pow25To7 = 6103515625.0; /* pow(25, 7) */
/*
* Step 1
*/
/* Equation 2 */
double C1 = sqrt((lab1.a * lab1.a) + (lab1.b * lab1.b));
double C2 = sqrt((lab2.a * lab2.a) + (lab2.b * lab2.b));
/* Equation 3 */
double barC = (C1 + C2) / 2.0;
/* Equation 4 */
double G = 0.5 * (1 - sqrt(pow(barC, 7) / (pow(barC, 7) + pow25To7)));
/* Equation 5 */
double a1Prime = (1.0 + G) * lab1.a;
double a2Prime = (1.0 + G) * lab2.a;
/* Equation 6 */
double CPrime1 = sqrt((a1Prime * a1Prime) + (lab1.b * lab1.b));
double CPrime2 = sqrt((a2Prime * a2Prime) + (lab2.b * lab2.b));
/* Equation 7 */
double hPrime1;
if (lab1.b == 0 && a1Prime == 0)
hPrime1 = 0.0;
else {
hPrime1 = atan2(lab1.b, a1Prime);
/*
* This must be converted to a hue angle in degrees between 0
* and 360 by addition of 2 to negative hue angles.
*/
if (hPrime1 < 0)
hPrime1 += deg360InRad;
}
double hPrime2;
if (lab2.b == 0 && a2Prime == 0)
hPrime2 = 0.0;
else {
hPrime2 = atan2(lab2.b, a2Prime);
/*
* This must be converted to a hue angle in degrees between 0
* and 360 by addition of 2 to negative hue angles.
*/
if (hPrime2 < 0)
hPrime2 += deg360InRad;
}
/*
* Step 2
*/
/* Equation 8 */
double deltaLPrime = lab2.l - lab1.l;
/* Equation 9 */
double deltaCPrime = CPrime2 - CPrime1;
/* Equation 10 */
double deltahPrime;
double CPrimeProduct = CPrime1 * CPrime2;
if (CPrimeProduct == 0)
deltahPrime = 0;
else {
/* Avoid the fabs() call */
deltahPrime = hPrime2 - hPrime1;
if (deltahPrime < -deg180InRad)
deltahPrime += deg360InRad;
else if (deltahPrime > deg180InRad)
deltahPrime -= deg360InRad;
}
/* Equation 11 */
double deltaHPrime = 2.0 * sqrt(CPrimeProduct) *
sin(deltahPrime / 2.0);
/*
* Step 3
*/
/* Equation 12 */
double barLPrime = (lab1.l + lab2.l) / 2.0;
/* Equation 13 */
double barCPrime = (CPrime1 + CPrime2) / 2.0;
/* Equation 14 */
double barhPrime, hPrimeSum = hPrime1 + hPrime2;
if (CPrime1 * CPrime2 == 0) {
barhPrime = hPrimeSum;
} else {
if (fabs(hPrime1 - hPrime2) <= deg180InRad)
barhPrime = hPrimeSum / 2.0;
else {
if (hPrimeSum < deg360InRad)
barhPrime = (hPrimeSum + deg360InRad) / 2.0;
else
barhPrime = (hPrimeSum - deg360InRad) / 2.0;
}
}
/* Equation 15 */
double T = 1.0 - (0.17 * cos(barhPrime - CIEDE2000::deg2Rad(30.0))) +
(0.24 * cos(2.0 * barhPrime)) +
(0.32 * cos((3.0 * barhPrime) + CIEDE2000::deg2Rad(6.0))) -
(0.20 * cos((4.0 * barhPrime) - CIEDE2000::deg2Rad(63.0)));
/* Equation 16 */
double deltaTheta = CIEDE2000::deg2Rad(30.0) *
exp(-pow((barhPrime - deg2Rad(275.0)) / deg2Rad(25.0), 2.0));
/* Equation 17 */
double R_C = 2.0 * sqrt(pow(barCPrime, 7.0) /
(pow(barCPrime, 7.0) + pow25To7));
/* Equation 18 */
double S_L = 1 + ((0.015 * pow(barLPrime - 50.0, 2.0)) /
sqrt(20 + pow(barLPrime - 50.0, 2.0)));
/* Equation 19 */
double S_C = 1 + (0.045 * barCPrime);
/* Equation 20 */
double S_H = 1 + (0.015 * barCPrime * T);
/* Equation 21 */
double R_T = (-sin(2.0 * deltaTheta)) * R_C;
/* Equation 22 */
double deltaE = sqrt(
pow(deltaLPrime / (k_L * S_L), 2.0) +
pow(deltaCPrime / (k_C * S_C), 2.0) +
pow(deltaHPrime / (k_H * S_H), 2.0) +
(R_T * (deltaCPrime / (k_C * S_C)) * (deltaHPrime / (k_H * S_H))));
return (deltaE);
}
/*******************************************************************************
* Operators.
******************************************************************************/
std::ostream&
operator<<(
std::ostream &s,
const CIEDE2000::LAB &labColor)
{
return (s << "CIELAB(" << labColor.l << "," << labColor.a << "," <<
labColor.b << ")");
}