1 //---------------------------------------------------------------------------------
2 //
3 //  Little Color Management System
4 //  Copyright (c) 1998-2010 Marti Maria Saguer
5 //
6 // Permission is hereby granted, free of charge, to any person obtaining
7 // a copy of this software and associated documentation files (the "Software"),
8 // to deal in the Software without restriction, including without limitation
9 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
10 // and/or sell copies of the Software, and to permit persons to whom the Software
11 // is furnished to do so, subject to the following conditions:
12 //
13 // The above copyright notice and this permission notice shall be included in
14 // all copies or substantial portions of the Software.
15 //
16 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
18 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
20 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
21 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
22 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23 //
24 //---------------------------------------------------------------------------------
25 //
26 
27 #include "lcms2_internal.h"
28 
29 //      inter PCS conversions XYZ <-> CIE L* a* b*
30 /*
31 
32 
33        CIE 15:2004 CIELab is defined as:
34 
35        L* = 116*f(Y/Yn) - 16                     0 <= L* <= 100
36        a* = 500*[f(X/Xn) - f(Y/Yn)]
37        b* = 200*[f(Y/Yn) - f(Z/Zn)]
38 
39        and
40 
41               f(t) = t^(1/3)                     1 >= t >  (24/116)^3
42                      (841/108)*t + (16/116)      0 <= t <= (24/116)^3
43 
44 
45        Reverse transform is:
46 
47        X = Xn*[a* / 500 + (L* + 16) / 116] ^ 3   if (X/Xn) > (24/116)
48          = Xn*(a* / 500 + L* / 116) / 7.787      if (X/Xn) <= (24/116)
49 
50 
51 
52        PCS in Lab2 is encoded as:
53 
54               8 bit Lab PCS:
55 
56                      L*      0..100 into a 0..ff byte.
57                      a*      t + 128 range is -128.0  +127.0
58                      b*
59 
60              16 bit Lab PCS:
61 
62                      L*     0..100  into a 0..ff00 word.
63                      a*     t + 128  range is  -128.0  +127.9961
64                      b*
65 
66 
67 
68 Interchange Space   Component     Actual Range        Encoded Range
69 CIE XYZ             X             0 -> 1.99997        0x0000 -> 0xffff
70 CIE XYZ             Y             0 -> 1.99997        0x0000 -> 0xffff
71 CIE XYZ             Z             0 -> 1.99997        0x0000 -> 0xffff
72 
73 Version 2,3
74 -----------
75 
76 CIELAB (16 bit)     L*            0 -> 100.0          0x0000 -> 0xff00
77 CIELAB (16 bit)     a*            -128.0 -> +127.996  0x0000 -> 0x8000 -> 0xffff
78 CIELAB (16 bit)     b*            -128.0 -> +127.996  0x0000 -> 0x8000 -> 0xffff
79 
80 
81 Version 4
82 ---------
83 
84 CIELAB (16 bit)     L*            0 -> 100.0          0x0000 -> 0xffff
85 CIELAB (16 bit)     a*            -128.0 -> +127      0x0000 -> 0x8080 -> 0xffff
86 CIELAB (16 bit)     b*            -128.0 -> +127      0x0000 -> 0x8080 -> 0xffff
87 
88 */
89 
90 // Conversions
cmsXYZ2xyY(cmsCIExyY * Dest,const cmsCIEXYZ * Source)91 void CMSEXPORT cmsXYZ2xyY(cmsCIExyY* Dest, const cmsCIEXYZ* Source)
92 {
93     cmsFloat64Number ISum;
94 
95     ISum = 1./(Source -> X + Source -> Y + Source -> Z);
96 
97     Dest -> x = (Source -> X) * ISum;
98     Dest -> y = (Source -> Y) * ISum;
99     Dest -> Y = Source -> Y;
100 }
101 
cmsxyY2XYZ(cmsCIEXYZ * Dest,const cmsCIExyY * Source)102 void CMSEXPORT cmsxyY2XYZ(cmsCIEXYZ* Dest, const cmsCIExyY* Source)
103 {
104     Dest -> X = (Source -> x / Source -> y) * Source -> Y;
105     Dest -> Y = Source -> Y;
106     Dest -> Z = ((1 - Source -> x - Source -> y) / Source -> y) * Source -> Y;
107 }
108 
109 static
f(cmsFloat64Number t)110 cmsFloat64Number f(cmsFloat64Number t)
111 {
112     const cmsFloat64Number Limit = (24.0/116.0) * (24.0/116.0) * (24.0/116.0);
113 
114     if (t <= Limit)
115         return (841.0/108.0) * t + (16.0/116.0);
116     else
117         return pow(t, 1.0/3.0);
118 }
119 
120 static
f_1(cmsFloat64Number t)121 cmsFloat64Number f_1(cmsFloat64Number t)
122 {
123     const cmsFloat64Number Limit = (24.0/116.0);
124 
125     if (t <= Limit) {
126         return (108.0/841.0) * (t - (16.0/116.0));
127     }
128 
129     return t * t * t;
130 }
131 
132 
133 // Standard XYZ to Lab. it can handle negative XZY numbers in some cases
cmsXYZ2Lab(const cmsCIEXYZ * WhitePoint,cmsCIELab * Lab,const cmsCIEXYZ * xyz)134 void CMSEXPORT cmsXYZ2Lab(const cmsCIEXYZ* WhitePoint, cmsCIELab* Lab, const cmsCIEXYZ* xyz)
135 {
136     cmsFloat64Number fx, fy, fz;
137 
138     if (WhitePoint == NULL)
139         WhitePoint = cmsD50_XYZ();
140 
141     fx = f(xyz->X / WhitePoint->X);
142     fy = f(xyz->Y / WhitePoint->Y);
143     fz = f(xyz->Z / WhitePoint->Z);
144 
145     Lab->L = 116.0*fy - 16.0;
146     Lab->a = 500.0*(fx - fy);
147     Lab->b = 200.0*(fy - fz);
148 }
149 
150 
151 // Standard XYZ to Lab. It can return negative XYZ in some cases
cmsLab2XYZ(const cmsCIEXYZ * WhitePoint,cmsCIEXYZ * xyz,const cmsCIELab * Lab)152 void CMSEXPORT cmsLab2XYZ(const cmsCIEXYZ* WhitePoint, cmsCIEXYZ* xyz,  const cmsCIELab* Lab)
153 {
154     cmsFloat64Number x, y, z;
155 
156     if (WhitePoint == NULL)
157         WhitePoint = cmsD50_XYZ();
158 
159     y = (Lab-> L + 16.0) / 116.0;
160     x = y + 0.002 * Lab -> a;
161     z = y - 0.005 * Lab -> b;
162 
163     xyz -> X = f_1(x) * WhitePoint -> X;
164     xyz -> Y = f_1(y) * WhitePoint -> Y;
165     xyz -> Z = f_1(z) * WhitePoint -> Z;
166 
167 }
168 
169 static
L2float2(cmsUInt16Number v)170 cmsFloat64Number L2float2(cmsUInt16Number v)
171 {
172     return (cmsFloat64Number) v / 652.800;
173 }
174 
175 // the a/b part
176 static
ab2float2(cmsUInt16Number v)177 cmsFloat64Number ab2float2(cmsUInt16Number v)
178 {
179     return ((cmsFloat64Number) v / 256.0) - 128.0;
180 }
181 
182 static
L2Fix2(cmsFloat64Number L)183 cmsUInt16Number L2Fix2(cmsFloat64Number L)
184 {
185     return _cmsQuickSaturateWord(L *  652.8);
186 }
187 
188 static
ab2Fix2(cmsFloat64Number ab)189 cmsUInt16Number ab2Fix2(cmsFloat64Number ab)
190 {
191     return _cmsQuickSaturateWord((ab + 128.0) * 256.0);
192 }
193 
194 
195 static
L2float4(cmsUInt16Number v)196 cmsFloat64Number L2float4(cmsUInt16Number v)
197 {
198     return (cmsFloat64Number) v / 655.35;
199 }
200 
201 // the a/b part
202 static
ab2float4(cmsUInt16Number v)203 cmsFloat64Number ab2float4(cmsUInt16Number v)
204 {
205     return ((cmsFloat64Number) v / 257.0) - 128.0;
206 }
207 
208 
cmsLabEncoded2FloatV2(cmsCIELab * Lab,const cmsUInt16Number wLab[3])209 void CMSEXPORT cmsLabEncoded2FloatV2(cmsCIELab* Lab, const cmsUInt16Number wLab[3])
210 {
211         Lab->L = L2float2(wLab[0]);
212         Lab->a = ab2float2(wLab[1]);
213         Lab->b = ab2float2(wLab[2]);
214 }
215 
216 
cmsLabEncoded2Float(cmsCIELab * Lab,const cmsUInt16Number wLab[3])217 void CMSEXPORT cmsLabEncoded2Float(cmsCIELab* Lab, const cmsUInt16Number wLab[3])
218 {
219         Lab->L = L2float4(wLab[0]);
220         Lab->a = ab2float4(wLab[1]);
221         Lab->b = ab2float4(wLab[2]);
222 }
223 
224 static
Clamp_L_doubleV2(cmsFloat64Number L)225 cmsFloat64Number Clamp_L_doubleV2(cmsFloat64Number L)
226 {
227     const cmsFloat64Number L_max = (cmsFloat64Number) (0xFFFF * 100.0) / 0xFF00;
228 
229     if (L < 0) L = 0;
230     if (L > L_max) L = L_max;
231 
232     return L;
233 }
234 
235 
236 static
Clamp_ab_doubleV2(cmsFloat64Number ab)237 cmsFloat64Number Clamp_ab_doubleV2(cmsFloat64Number ab)
238 {
239     if (ab < MIN_ENCODEABLE_ab2) ab = MIN_ENCODEABLE_ab2;
240     if (ab > MAX_ENCODEABLE_ab2) ab = MAX_ENCODEABLE_ab2;
241 
242     return ab;
243 }
244 
cmsFloat2LabEncodedV2(cmsUInt16Number wLab[3],const cmsCIELab * fLab)245 void CMSEXPORT cmsFloat2LabEncodedV2(cmsUInt16Number wLab[3], const cmsCIELab* fLab)
246 {
247     cmsCIELab Lab;
248 
249     Lab.L = Clamp_L_doubleV2(fLab ->L);
250     Lab.a = Clamp_ab_doubleV2(fLab ->a);
251     Lab.b = Clamp_ab_doubleV2(fLab ->b);
252 
253     wLab[0] = L2Fix2(Lab.L);
254     wLab[1] = ab2Fix2(Lab.a);
255     wLab[2] = ab2Fix2(Lab.b);
256 }
257 
258 
259 static
Clamp_L_doubleV4(cmsFloat64Number L)260 cmsFloat64Number Clamp_L_doubleV4(cmsFloat64Number L)
261 {
262     if (L < 0) L = 0;
263     if (L > 100.0) L = 100.0;
264 
265     return L;
266 }
267 
268 static
Clamp_ab_doubleV4(cmsFloat64Number ab)269 cmsFloat64Number Clamp_ab_doubleV4(cmsFloat64Number ab)
270 {
271     if (ab < MIN_ENCODEABLE_ab4) ab = MIN_ENCODEABLE_ab4;
272     if (ab > MAX_ENCODEABLE_ab4) ab = MAX_ENCODEABLE_ab4;
273 
274     return ab;
275 }
276 
277 static
L2Fix4(cmsFloat64Number L)278 cmsUInt16Number L2Fix4(cmsFloat64Number L)
279 {
280     return _cmsQuickSaturateWord(L *  655.35);
281 }
282 
283 static
ab2Fix4(cmsFloat64Number ab)284 cmsUInt16Number ab2Fix4(cmsFloat64Number ab)
285 {
286     return _cmsQuickSaturateWord((ab + 128.0) * 257.0);
287 }
288 
cmsFloat2LabEncoded(cmsUInt16Number wLab[3],const cmsCIELab * fLab)289 void CMSEXPORT cmsFloat2LabEncoded(cmsUInt16Number wLab[3], const cmsCIELab* fLab)
290 {
291     cmsCIELab Lab;
292 
293     Lab.L = Clamp_L_doubleV4(fLab ->L);
294     Lab.a = Clamp_ab_doubleV4(fLab ->a);
295     Lab.b = Clamp_ab_doubleV4(fLab ->b);
296 
297     wLab[0] = L2Fix4(Lab.L);
298     wLab[1] = ab2Fix4(Lab.a);
299     wLab[2] = ab2Fix4(Lab.b);
300 }
301 
302 // Auxiliar: convert to Radians
303 static
RADIANS(cmsFloat64Number deg)304 cmsFloat64Number RADIANS(cmsFloat64Number deg)
305 {
306     return (deg * M_PI) / 180.;
307 }
308 
309 
310 // Auxiliar: atan2 but operating in degrees and returning 0 if a==b==0
311 static
atan2deg(cmsFloat64Number a,cmsFloat64Number b)312 cmsFloat64Number atan2deg(cmsFloat64Number a, cmsFloat64Number b)
313 {
314    cmsFloat64Number h;
315 
316    if (a == 0 && b == 0)
317             h   = 0;
318     else
319             h = atan2(a, b);
320 
321     h *= (180. / M_PI);
322 
323     while (h > 360.)
324         h -= 360.;
325 
326     while ( h < 0)
327         h += 360.;
328 
329     return h;
330 }
331 
332 
333 // Auxiliar: Square
334 static
Sqr(cmsFloat64Number v)335 cmsFloat64Number Sqr(cmsFloat64Number v)
336 {
337     return v *  v;
338 }
339 // From cylindrical coordinates. No check is performed, then negative values are allowed
cmsLab2LCh(cmsCIELCh * LCh,const cmsCIELab * Lab)340 void CMSEXPORT cmsLab2LCh(cmsCIELCh* LCh, const cmsCIELab* Lab)
341 {
342     LCh -> L = Lab -> L;
343     LCh -> C = pow(Sqr(Lab ->a) + Sqr(Lab ->b), 0.5);
344     LCh -> h = atan2deg(Lab ->b, Lab ->a);
345 }
346 
347 
348 // To cylindrical coordinates. No check is performed, then negative values are allowed
cmsLCh2Lab(cmsCIELab * Lab,const cmsCIELCh * LCh)349 void CMSEXPORT cmsLCh2Lab(cmsCIELab* Lab, const cmsCIELCh* LCh)
350 {
351     cmsFloat64Number h = (LCh -> h * M_PI) / 180.0;
352 
353     Lab -> L = LCh -> L;
354     Lab -> a = LCh -> C * cos(h);
355     Lab -> b = LCh -> C * sin(h);
356 }
357 
358 // In XYZ All 3 components are encoded using 1.15 fixed point
359 static
XYZ2Fix(cmsFloat64Number d)360 cmsUInt16Number XYZ2Fix(cmsFloat64Number d)
361 {
362     return _cmsQuickSaturateWord(d * 32768.0);
363 }
364 
cmsFloat2XYZEncoded(cmsUInt16Number XYZ[3],const cmsCIEXYZ * fXYZ)365 void CMSEXPORT cmsFloat2XYZEncoded(cmsUInt16Number XYZ[3], const cmsCIEXYZ* fXYZ)
366 {
367     cmsCIEXYZ xyz;
368 
369     xyz.X = fXYZ -> X;
370     xyz.Y = fXYZ -> Y;
371     xyz.Z = fXYZ -> Z;
372 
373     // Clamp to encodeable values.
374     if (xyz.Y <= 0) {
375 
376         xyz.X = 0;
377         xyz.Y = 0;
378         xyz.Z = 0;
379     }
380 
381     if (xyz.X > MAX_ENCODEABLE_XYZ)
382         xyz.X = MAX_ENCODEABLE_XYZ;
383 
384     if (xyz.X < 0)
385         xyz.X = 0;
386 
387     if (xyz.Y > MAX_ENCODEABLE_XYZ)
388         xyz.Y = MAX_ENCODEABLE_XYZ;
389 
390     if (xyz.Y < 0)
391         xyz.Y = 0;
392 
393     if (xyz.Z > MAX_ENCODEABLE_XYZ)
394         xyz.Z = MAX_ENCODEABLE_XYZ;
395 
396     if (xyz.Z < 0)
397         xyz.Z = 0;
398 
399 
400     XYZ[0] = XYZ2Fix(xyz.X);
401     XYZ[1] = XYZ2Fix(xyz.Y);
402     XYZ[2] = XYZ2Fix(xyz.Z);
403 }
404 
405 
406 //  To convert from Fixed 1.15 point to cmsFloat64Number
407 static
XYZ2float(cmsUInt16Number v)408 cmsFloat64Number XYZ2float(cmsUInt16Number v)
409 {
410     cmsS15Fixed16Number fix32;
411 
412     // From 1.15 to 15.16
413     fix32 = v << 1;
414 
415     // From fixed 15.16 to cmsFloat64Number
416     return _cms15Fixed16toDouble(fix32);
417 }
418 
419 
cmsXYZEncoded2Float(cmsCIEXYZ * fXYZ,const cmsUInt16Number XYZ[3])420 void CMSEXPORT cmsXYZEncoded2Float(cmsCIEXYZ* fXYZ, const cmsUInt16Number XYZ[3])
421 {
422     fXYZ -> X = XYZ2float(XYZ[0]);
423     fXYZ -> Y = XYZ2float(XYZ[1]);
424     fXYZ -> Z = XYZ2float(XYZ[2]);
425 }
426 
427 
428 // Returns dE on two Lab values
cmsDeltaE(const cmsCIELab * Lab1,const cmsCIELab * Lab2)429 cmsFloat64Number CMSEXPORT cmsDeltaE(const cmsCIELab* Lab1, const cmsCIELab* Lab2)
430 {
431     cmsFloat64Number dL, da, db;
432 
433     dL = fabs(Lab1 -> L - Lab2 -> L);
434     da = fabs(Lab1 -> a - Lab2 -> a);
435     db = fabs(Lab1 -> b - Lab2 -> b);
436 
437     return pow(Sqr(dL) + Sqr(da) + Sqr(db), 0.5);
438 }
439 
440 
441 // Return the CIE94 Delta E
cmsCIE94DeltaE(const cmsCIELab * Lab1,const cmsCIELab * Lab2)442 cmsFloat64Number CMSEXPORT cmsCIE94DeltaE(const cmsCIELab* Lab1, const cmsCIELab* Lab2)
443 {
444     cmsCIELCh LCh1, LCh2;
445     cmsFloat64Number dE, dL, dC, dh, dhsq;
446     cmsFloat64Number c12, sc, sh;
447 
448     dL = fabs(Lab1 ->L - Lab2 ->L);
449 
450     cmsLab2LCh(&LCh1, Lab1);
451     cmsLab2LCh(&LCh2, Lab2);
452 
453     dC  = fabs(LCh1.C - LCh2.C);
454     dE  = cmsDeltaE(Lab1, Lab2);
455 
456     dhsq = Sqr(dE) - Sqr(dL) - Sqr(dC);
457     if (dhsq < 0)
458         dh = 0;
459     else
460         dh = pow(dhsq, 0.5);
461 
462     c12 = sqrt(LCh1.C * LCh2.C);
463 
464     sc = 1.0 + (0.048 * c12);
465     sh = 1.0 + (0.014 * c12);
466 
467     return sqrt(Sqr(dL)  + Sqr(dC) / Sqr(sc) + Sqr(dh) / Sqr(sh));
468 }
469 
470 
471 // Auxiliary
472 static
ComputeLBFD(const cmsCIELab * Lab)473 cmsFloat64Number ComputeLBFD(const cmsCIELab* Lab)
474 {
475   cmsFloat64Number yt;
476 
477   if (Lab->L > 7.996969)
478         yt = (Sqr((Lab->L+16)/116)*((Lab->L+16)/116))*100;
479   else
480         yt = 100 * (Lab->L / 903.3);
481 
482   return (54.6 * (M_LOG10E * (log(yt + 1.5))) - 9.6);
483 }
484 
485 
486 
487 // bfd - gets BFD(1:1) difference between Lab1, Lab2
cmsBFDdeltaE(const cmsCIELab * Lab1,const cmsCIELab * Lab2)488 cmsFloat64Number CMSEXPORT cmsBFDdeltaE(const cmsCIELab* Lab1, const cmsCIELab* Lab2)
489 {
490     cmsFloat64Number lbfd1,lbfd2,AveC,Aveh,dE,deltaL,
491         deltaC,deltah,dc,t,g,dh,rh,rc,rt,bfd;
492     cmsCIELCh LCh1, LCh2;
493 
494 
495     lbfd1 = ComputeLBFD(Lab1);
496     lbfd2 = ComputeLBFD(Lab2);
497     deltaL = lbfd2 - lbfd1;
498 
499     cmsLab2LCh(&LCh1, Lab1);
500     cmsLab2LCh(&LCh2, Lab2);
501 
502     deltaC = LCh2.C - LCh1.C;
503     AveC = (LCh1.C+LCh2.C)/2;
504     Aveh = (LCh1.h+LCh2.h)/2;
505 
506     dE = cmsDeltaE(Lab1, Lab2);
507 
508     if (Sqr(dE)>(Sqr(Lab2->L-Lab1->L)+Sqr(deltaC)))
509         deltah = sqrt(Sqr(dE)-Sqr(Lab2->L-Lab1->L)-Sqr(deltaC));
510     else
511         deltah =0;
512 
513 
514     dc   = 0.035 * AveC / (1 + 0.00365 * AveC)+0.521;
515     g    = sqrt(Sqr(Sqr(AveC))/(Sqr(Sqr(AveC))+14000));
516     t    = 0.627+(0.055*cos((Aveh-254)/(180/M_PI))-
517            0.040*cos((2*Aveh-136)/(180/M_PI))+
518            0.070*cos((3*Aveh-31)/(180/M_PI))+
519            0.049*cos((4*Aveh+114)/(180/M_PI))-
520            0.015*cos((5*Aveh-103)/(180/M_PI)));
521 
522     dh    = dc*(g*t+1-g);
523     rh    = -0.260*cos((Aveh-308)/(180/M_PI))-
524            0.379*cos((2*Aveh-160)/(180/M_PI))-
525            0.636*cos((3*Aveh+254)/(180/M_PI))+
526            0.226*cos((4*Aveh+140)/(180/M_PI))-
527            0.194*cos((5*Aveh+280)/(180/M_PI));
528 
529     rc = sqrt((AveC*AveC*AveC*AveC*AveC*AveC)/((AveC*AveC*AveC*AveC*AveC*AveC)+70000000));
530     rt = rh*rc;
531 
532     bfd = sqrt(Sqr(deltaL)+Sqr(deltaC/dc)+Sqr(deltah/dh)+(rt*(deltaC/dc)*(deltah/dh)));
533 
534     return bfd;
535 }
536 
537 
538 //  cmc - CMC(l:c) difference between Lab1, Lab2
cmsCMCdeltaE(const cmsCIELab * Lab1,const cmsCIELab * Lab2,cmsFloat64Number l,cmsFloat64Number c)539 cmsFloat64Number CMSEXPORT cmsCMCdeltaE(const cmsCIELab* Lab1, const cmsCIELab* Lab2, cmsFloat64Number l, cmsFloat64Number c)
540 {
541   cmsFloat64Number dE,dL,dC,dh,sl,sc,sh,t,f,cmc;
542   cmsCIELCh LCh1, LCh2;
543 
544   if (Lab1 ->L == 0 && Lab2 ->L == 0) return 0;
545 
546   cmsLab2LCh(&LCh1, Lab1);
547   cmsLab2LCh(&LCh2, Lab2);
548 
549 
550   dL = Lab2->L-Lab1->L;
551   dC = LCh2.C-LCh1.C;
552 
553   dE = cmsDeltaE(Lab1, Lab2);
554 
555   if (Sqr(dE)>(Sqr(dL)+Sqr(dC)))
556             dh = sqrt(Sqr(dE)-Sqr(dL)-Sqr(dC));
557   else
558             dh =0;
559 
560   if ((LCh1.h > 164) && (LCh1.h < 345))
561       t = 0.56 + fabs(0.2 * cos(((LCh1.h + 168)/(180/M_PI))));
562   else
563       t = 0.36 + fabs(0.4 * cos(((LCh1.h + 35 )/(180/M_PI))));
564 
565    sc  = 0.0638   * LCh1.C / (1 + 0.0131  * LCh1.C) + 0.638;
566    sl  = 0.040975 * Lab1->L /(1 + 0.01765 * Lab1->L);
567 
568    if (Lab1->L<16)
569          sl = 0.511;
570 
571    f   = sqrt((LCh1.C * LCh1.C * LCh1.C * LCh1.C)/((LCh1.C * LCh1.C * LCh1.C * LCh1.C)+1900));
572    sh  = sc*(t*f+1-f);
573    cmc = sqrt(Sqr(dL/(l*sl))+Sqr(dC/(c*sc))+Sqr(dh/sh));
574 
575    return cmc;
576 }
577 
578 // dE2000 The weightings KL, KC and KH can be modified to reflect the relative
579 // importance of lightness, chroma and hue in different industrial applications
cmsCIE2000DeltaE(const cmsCIELab * Lab1,const cmsCIELab * Lab2,cmsFloat64Number Kl,cmsFloat64Number Kc,cmsFloat64Number Kh)580 cmsFloat64Number CMSEXPORT cmsCIE2000DeltaE(const cmsCIELab* Lab1, const cmsCIELab* Lab2,
581                                   cmsFloat64Number Kl, cmsFloat64Number Kc, cmsFloat64Number Kh)
582 {
583     cmsFloat64Number L1  = Lab1->L;
584     cmsFloat64Number a1  = Lab1->a;
585     cmsFloat64Number b1  = Lab1->b;
586     cmsFloat64Number C   = sqrt( Sqr(a1) + Sqr(b1) );
587 
588     cmsFloat64Number Ls = Lab2 ->L;
589     cmsFloat64Number as = Lab2 ->a;
590     cmsFloat64Number bs = Lab2 ->b;
591     cmsFloat64Number Cs = sqrt( Sqr(as) + Sqr(bs) );
592 
593     cmsFloat64Number G = 0.5 * ( 1 - sqrt(pow((C + Cs) / 2 , 7.0) / (pow((C + Cs) / 2, 7.0) + pow(25.0, 7.0) ) ));
594 
595     cmsFloat64Number a_p = (1 + G ) * a1;
596     cmsFloat64Number b_p = b1;
597     cmsFloat64Number C_p = sqrt( Sqr(a_p) + Sqr(b_p));
598     cmsFloat64Number h_p = atan2deg(b_p, a_p);
599 
600 
601     cmsFloat64Number a_ps = (1 + G) * as;
602     cmsFloat64Number b_ps = bs;
603     cmsFloat64Number C_ps = sqrt(Sqr(a_ps) + Sqr(b_ps));
604     cmsFloat64Number h_ps = atan2deg(b_ps, a_ps);
605 
606     cmsFloat64Number meanC_p =(C_p + C_ps) / 2;
607 
608     cmsFloat64Number hps_plus_hp  = h_ps + h_p;
609     cmsFloat64Number hps_minus_hp = h_ps - h_p;
610 
611     cmsFloat64Number meanh_p = fabs(hps_minus_hp) <= 180.000001 ? (hps_plus_hp)/2 :
612                             (hps_plus_hp) < 360 ? (hps_plus_hp + 360)/2 :
613                                                  (hps_plus_hp - 360)/2;
614 
615     cmsFloat64Number delta_h = (hps_minus_hp) <= -180.000001 ?  (hps_minus_hp + 360) :
616                             (hps_minus_hp) > 180 ? (hps_minus_hp - 360) :
617                                                     (hps_minus_hp);
618     cmsFloat64Number delta_L = (Ls - L1);
619     cmsFloat64Number delta_C = (C_ps - C_p );
620 
621 
622     cmsFloat64Number delta_H =2 * sqrt(C_ps*C_p) * sin(RADIANS(delta_h) / 2);
623 
624     cmsFloat64Number T = 1 - 0.17 * cos(RADIANS(meanh_p-30))
625                  + 0.24 * cos(RADIANS(2*meanh_p))
626                  + 0.32 * cos(RADIANS(3*meanh_p + 6))
627                  - 0.2  * cos(RADIANS(4*meanh_p - 63));
628 
629     cmsFloat64Number Sl = 1 + (0.015 * Sqr((Ls + L1) /2- 50) )/ sqrt(20 + Sqr( (Ls+L1)/2 - 50) );
630 
631     cmsFloat64Number Sc = 1 + 0.045 * (C_p + C_ps)/2;
632     cmsFloat64Number Sh = 1 + 0.015 * ((C_ps + C_p)/2) * T;
633 
634     cmsFloat64Number delta_ro = 30 * exp( -Sqr(((meanh_p - 275 ) / 25)));
635 
636     cmsFloat64Number Rc = 2 * sqrt(( pow(meanC_p, 7.0) )/( pow(meanC_p, 7.0) + pow(25.0, 7.0)));
637 
638     cmsFloat64Number Rt = -sin(2 * RADIANS(delta_ro)) * Rc;
639 
640     cmsFloat64Number deltaE00 = sqrt( Sqr(delta_L /(Sl * Kl)) +
641                             Sqr(delta_C/(Sc * Kc))  +
642                             Sqr(delta_H/(Sh * Kh))  +
643                             Rt*(delta_C/(Sc * Kc)) * (delta_H / (Sh * Kh)));
644 
645     return deltaE00;
646 }
647 
648 // This function returns a number of gridpoints to be used as LUT table. It assumes same number
649 // of gripdpoints in all dimensions. Flags may override the choice.
_cmsReasonableGridpointsByColorspace(cmsColorSpaceSignature Colorspace,cmsUInt32Number dwFlags)650 int _cmsReasonableGridpointsByColorspace(cmsColorSpaceSignature Colorspace, cmsUInt32Number dwFlags)
651 {
652     int nChannels;
653 
654     // Already specified?
655     if (dwFlags & 0x00FF0000) {
656             // Yes, grab'em
657             return (dwFlags >> 16) & 0xFF;
658     }
659 
660     nChannels = cmsChannelsOf(Colorspace);
661 
662     // HighResPrecalc is maximum resolution
663     if (dwFlags & cmsFLAGS_HIGHRESPRECALC) {
664 
665         if (nChannels > 4)
666                 return 7;       // 7 for Hifi
667 
668         if (nChannels == 4)     // 23 for CMYK
669                 return 23;
670 
671         return 49;      // 49 for RGB and others
672     }
673 
674 
675     // LowResPrecal is lower resolution
676     if (dwFlags & cmsFLAGS_LOWRESPRECALC) {
677 
678         if (nChannels > 4)
679                 return 6;       // 6 for more than 4 channels
680 
681         if (nChannels == 1)
682                 return 33;      // For monochrome
683 
684         return 17;              // 17 for remaining
685     }
686 
687     // Default values
688     if (nChannels > 4)
689                 return 7;       // 7 for Hifi
690 
691     if (nChannels == 4)
692                 return 17;      // 17 for CMYK
693 
694     return 33;                  // 33 for RGB
695 }
696 
697 
_cmsEndPointsBySpace(cmsColorSpaceSignature Space,cmsUInt16Number ** White,cmsUInt16Number ** Black,cmsUInt32Number * nOutputs)698 cmsBool  _cmsEndPointsBySpace(cmsColorSpaceSignature Space,
699                              cmsUInt16Number **White,
700                              cmsUInt16Number **Black,
701                              cmsUInt32Number *nOutputs)
702 {
703        // Only most common spaces
704 
705        static cmsUInt16Number RGBblack[4]  = { 0, 0, 0 };
706        static cmsUInt16Number RGBwhite[4]  = { 0xffff, 0xffff, 0xffff };
707        static cmsUInt16Number CMYKblack[4] = { 0xffff, 0xffff, 0xffff, 0xffff };   // 400% of ink
708        static cmsUInt16Number CMYKwhite[4] = { 0, 0, 0, 0 };
709        static cmsUInt16Number LABblack[4]  = { 0, 0x8080, 0x8080 };               // V4 Lab encoding
710        static cmsUInt16Number LABwhite[4]  = { 0xFFFF, 0x8080, 0x8080 };
711        static cmsUInt16Number CMYblack[4]  = { 0xffff, 0xffff, 0xffff };
712        static cmsUInt16Number CMYwhite[4]  = { 0, 0, 0 };
713        static cmsUInt16Number Grayblack[4] = { 0 };
714        static cmsUInt16Number GrayWhite[4] = { 0xffff };
715 
716        switch (Space) {
717 
718        case cmsSigGrayData: if (White)    *White = GrayWhite;
719                            if (Black)    *Black = Grayblack;
720                            if (nOutputs) *nOutputs = 1;
721                            return TRUE;
722 
723        case cmsSigRgbData:  if (White)    *White = RGBwhite;
724                            if (Black)    *Black = RGBblack;
725                            if (nOutputs) *nOutputs = 3;
726                            return TRUE;
727 
728        case cmsSigLabData:  if (White)    *White = LABwhite;
729                            if (Black)    *Black = LABblack;
730                            if (nOutputs) *nOutputs = 3;
731                            return TRUE;
732 
733        case cmsSigCmykData: if (White)    *White = CMYKwhite;
734                            if (Black)    *Black = CMYKblack;
735                            if (nOutputs) *nOutputs = 4;
736                            return TRUE;
737 
738        case cmsSigCmyData:  if (White)    *White = CMYwhite;
739                            if (Black)    *Black = CMYblack;
740                            if (nOutputs) *nOutputs = 3;
741                            return TRUE;
742 
743        default:;
744        }
745 
746   return FALSE;
747 }
748 
749 
750 
751 // Several utilities -------------------------------------------------------
752 
753 // Translate from our colorspace to ICC representation
754 
_cmsICCcolorSpace(int OurNotation)755 cmsColorSpaceSignature CMSEXPORT _cmsICCcolorSpace(int OurNotation)
756 {
757        switch (OurNotation) {
758 
759        case 1:
760        case PT_GRAY: return cmsSigGrayData;
761 
762        case 2:
763        case PT_RGB:  return cmsSigRgbData;
764 
765        case PT_CMY:  return cmsSigCmyData;
766        case PT_CMYK: return cmsSigCmykData;
767        case PT_YCbCr:return cmsSigYCbCrData;
768        case PT_YUV:  return cmsSigLuvData;
769        case PT_XYZ:  return cmsSigXYZData;
770 
771        case PT_LabV2:
772        case PT_Lab:  return cmsSigLabData;
773 
774        case PT_YUVK: return cmsSigLuvKData;
775        case PT_HSV:  return cmsSigHsvData;
776        case PT_HLS:  return cmsSigHlsData;
777        case PT_Yxy:  return cmsSigYxyData;
778 
779        case PT_MCH1: return cmsSigMCH1Data;
780        case PT_MCH2: return cmsSigMCH2Data;
781        case PT_MCH3: return cmsSigMCH3Data;
782        case PT_MCH4: return cmsSigMCH4Data;
783        case PT_MCH5: return cmsSigMCH5Data;
784        case PT_MCH6: return cmsSigMCH6Data;
785        case PT_MCH7: return cmsSigMCH7Data;
786        case PT_MCH8: return cmsSigMCH8Data;
787 
788        case PT_MCH9:  return cmsSigMCH9Data;
789        case PT_MCH10: return cmsSigMCHAData;
790        case PT_MCH11: return cmsSigMCHBData;
791        case PT_MCH12: return cmsSigMCHCData;
792        case PT_MCH13: return cmsSigMCHDData;
793        case PT_MCH14: return cmsSigMCHEData;
794        case PT_MCH15: return cmsSigMCHFData;
795 
796        default:  return (cmsColorSpaceSignature) (-1);
797        }
798 }
799 
800 
_cmsLCMScolorSpace(cmsColorSpaceSignature ProfileSpace)801 int CMSEXPORT _cmsLCMScolorSpace(cmsColorSpaceSignature ProfileSpace)
802 {
803     switch (ProfileSpace) {
804 
805     case cmsSigGrayData: return  PT_GRAY;
806     case cmsSigRgbData:  return  PT_RGB;
807     case cmsSigCmyData:  return  PT_CMY;
808     case cmsSigCmykData: return  PT_CMYK;
809     case cmsSigYCbCrData:return  PT_YCbCr;
810     case cmsSigLuvData:  return  PT_YUV;
811     case cmsSigXYZData:  return  PT_XYZ;
812     case cmsSigLabData:  return  PT_Lab;
813     case cmsSigLuvKData: return  PT_YUVK;
814     case cmsSigHsvData:  return  PT_HSV;
815     case cmsSigHlsData:  return  PT_HLS;
816     case cmsSigYxyData:  return  PT_Yxy;
817 
818     case cmsSig1colorData:
819     case cmsSigMCH1Data: return PT_MCH1;
820 
821     case cmsSig2colorData:
822     case cmsSigMCH2Data: return PT_MCH2;
823 
824     case cmsSig3colorData:
825     case cmsSigMCH3Data: return PT_MCH3;
826 
827     case cmsSig4colorData:
828     case cmsSigMCH4Data: return PT_MCH4;
829 
830     case cmsSig5colorData:
831     case cmsSigMCH5Data: return PT_MCH5;
832 
833     case cmsSig6colorData:
834     case cmsSigMCH6Data: return PT_MCH6;
835 
836     case cmsSigMCH7Data:
837     case cmsSig7colorData:return PT_MCH7;
838 
839     case cmsSigMCH8Data:
840     case cmsSig8colorData:return PT_MCH8;
841 
842     case cmsSigMCH9Data:
843     case cmsSig9colorData:return PT_MCH9;
844 
845     case cmsSigMCHAData:
846     case cmsSig10colorData:return PT_MCH10;
847 
848     case cmsSigMCHBData:
849     case cmsSig11colorData:return PT_MCH11;
850 
851     case cmsSigMCHCData:
852     case cmsSig12colorData:return PT_MCH12;
853 
854     case cmsSigMCHDData:
855     case cmsSig13colorData:return PT_MCH13;
856 
857     case cmsSigMCHEData:
858     case cmsSig14colorData:return PT_MCH14;
859 
860     case cmsSigMCHFData:
861     case cmsSig15colorData:return PT_MCH15;
862 
863     default:  return (cmsColorSpaceSignature) (-1);
864     }
865 }
866 
867 
cmsChannelsOf(cmsColorSpaceSignature ColorSpace)868 cmsUInt32Number CMSEXPORT cmsChannelsOf(cmsColorSpaceSignature ColorSpace)
869 {
870     switch (ColorSpace) {
871 
872     case cmsSigMCH1Data:
873     case cmsSig1colorData:
874     case cmsSigGrayData: return 1;
875 
876     case cmsSigMCH2Data:
877     case cmsSig2colorData:  return 2;
878 
879     case cmsSigXYZData:
880     case cmsSigLabData:
881     case cmsSigLuvData:
882     case cmsSigYCbCrData:
883     case cmsSigYxyData:
884     case cmsSigRgbData:
885     case cmsSigHsvData:
886     case cmsSigHlsData:
887     case cmsSigCmyData:
888     case cmsSigMCH3Data:
889     case cmsSig3colorData:  return 3;
890 
891     case cmsSigLuvKData:
892     case cmsSigCmykData:
893     case cmsSigMCH4Data:
894     case cmsSig4colorData:  return 4;
895 
896     case cmsSigMCH5Data:
897     case cmsSig5colorData:  return 5;
898 
899     case cmsSigMCH6Data:
900     case cmsSig6colorData:  return 6;
901 
902     case cmsSigMCH7Data:
903     case cmsSig7colorData:  return  7;
904 
905     case cmsSigMCH8Data:
906     case cmsSig8colorData:  return  8;
907 
908     case cmsSigMCH9Data:
909     case cmsSig9colorData:  return  9;
910 
911     case cmsSigMCHAData:
912     case cmsSig10colorData: return 10;
913 
914     case cmsSigMCHBData:
915     case cmsSig11colorData: return 11;
916 
917     case cmsSigMCHCData:
918     case cmsSig12colorData: return 12;
919 
920     case cmsSigMCHDData:
921     case cmsSig13colorData: return 13;
922 
923     case cmsSigMCHEData:
924     case cmsSig14colorData: return 14;
925 
926     case cmsSigMCHFData:
927     case cmsSig15colorData: return 15;
928 
929     default: return 3;
930     }
931 }
932