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
68Interchange Space   Component     Actual Range        Encoded Range
69CIE XYZ             X             0 -> 1.99997        0x0000 -> 0xffff
70CIE XYZ             Y             0 -> 1.99997        0x0000 -> 0xffff
71CIE XYZ             Z             0 -> 1.99997        0x0000 -> 0xffff
72
73Version 2,3
74-----------
75
76CIELAB (16 bit)     L*            0 -> 100.0          0x0000 -> 0xff00
77CIELAB (16 bit)     a*            -128.0 -> +127.996  0x0000 -> 0x8000 -> 0xffff
78CIELAB (16 bit)     b*            -128.0 -> +127.996  0x0000 -> 0x8000 -> 0xffff
79
80
81Version 4
82---------
83
84CIELAB (16 bit)     L*            0 -> 100.0          0x0000 -> 0xffff
85CIELAB (16 bit)     a*            -128.0 -> +127      0x0000 -> 0x8080 -> 0xffff
86CIELAB (16 bit)     b*            -128.0 -> +127      0x0000 -> 0x8080 -> 0xffff
87
88*/
89
90// Conversions
91void 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
102void 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
109static
110cmsFloat64Number 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
120static
121cmsFloat64Number 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
134void 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
152void 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
169static
170cmsFloat64Number L2float2(cmsUInt16Number v)
171{
172    return (cmsFloat64Number) v / 652.800;
173}
174
175// the a/b part
176static
177cmsFloat64Number ab2float2(cmsUInt16Number v)
178{
179    return ((cmsFloat64Number) v / 256.0) - 128.0;
180}
181
182static
183cmsUInt16Number L2Fix2(cmsFloat64Number L)
184{
185    return _cmsQuickSaturateWord(L *  652.8);
186}
187
188static
189cmsUInt16Number ab2Fix2(cmsFloat64Number ab)
190{
191    return _cmsQuickSaturateWord((ab + 128.0) * 256.0);
192}
193
194
195static
196cmsFloat64Number L2float4(cmsUInt16Number v)
197{
198    return (cmsFloat64Number) v / 655.35;
199}
200
201// the a/b part
202static
203cmsFloat64Number ab2float4(cmsUInt16Number v)
204{
205    return ((cmsFloat64Number) v / 257.0) - 128.0;
206}
207
208
209void 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
217void 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
224static
225cmsFloat64Number 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
236static
237cmsFloat64Number 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
245void 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
259static
260cmsFloat64Number 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
268static
269cmsFloat64Number 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
277static
278cmsUInt16Number L2Fix4(cmsFloat64Number L)
279{
280    return _cmsQuickSaturateWord(L *  655.35);
281}
282
283static
284cmsUInt16Number ab2Fix4(cmsFloat64Number ab)
285{
286    return _cmsQuickSaturateWord((ab + 128.0) * 257.0);
287}
288
289void 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
303static
304cmsFloat64Number 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
311static
312cmsFloat64Number 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
334static
335cmsFloat64Number Sqr(cmsFloat64Number v)
336{
337    return v *  v;
338}
339// From cylindrical coordinates. No check is performed, then negative values are allowed
340void 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
349void 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
359static
360cmsUInt16Number XYZ2Fix(cmsFloat64Number d)
361{
362    return _cmsQuickSaturateWord(d * 32768.0);
363}
364
365void 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
407static
408cmsFloat64Number 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
420void 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
429cmsFloat64Number 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
442cmsFloat64Number 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
472static
473cmsFloat64Number 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
488cmsFloat64Number 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
539cmsFloat64Number 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
580cmsFloat64Number 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.
650int _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
698cmsBool  _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
755cmsColorSpaceSignature 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
801int 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
868cmsUInt32Number 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