thirdparty/liblcms2/src/cmspcs.c

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