modules/lcms/src/cmspcs.c

   1 //
   2 //  Little cms
   3 //  Copyright (C) 1998-2007 Marti Maria
   4 //
   5 // Permission is hereby granted, free of charge, to any person obtaining
   6 // a copy of this software and associated documentation files (the "Software"),
   7 // to deal in the Software without restriction, including without limitation
   8 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
   9 // and/or sell copies of the Software, and to permit persons to whom the Software
  10 // is furnished to do so, subject to the following conditions:
  11 //
  12 // The above copyright notice and this permission notice shall be included in
  13 // all copies or substantial portions of the Software.
  14 //
  15 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  16 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
  17 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  18 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  19 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  20 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  21 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  22
  23 //      inter PCS conversions XYZ <-> CIE L* a* b*
  24
  25 #include "lcms.h"
  26
  27 /*
  28
  29
  30        CIE 15:2004 CIELab is defined as:
  31
  32        L* = 116*f(Y/Yn) - 16                     0 <= L* <= 100
  33        a* = 500*[f(X/Xn) - f(Y/Yn)]
  34        b* = 200*[f(Y/Yn) - f(Z/Zn)]
  35
  36        and
  37
  38               f(t) = t^(1/3)                     1 >= t >  (24/116)^3
  39                      (841/108)*t + (16/116)      0 <= t <= (24/116)^3
  40
  41
  42        Reverse transform is:
  43
  44        X = Xn*[a* / 500 + (L* + 16) / 116] ^ 3   if (X/Xn) > (24/116)
  45          = Xn*(a* / 500 + L* / 116) / 7.787      if (X/Xn) <= (24/116)
  46
  47
  48
  49        Following ICC. PCS in Lab is coded as:
  50
  51               8 bit Lab PCS:
  52
  53                      L*      0..100 into a 0..ff byte.
  54                      a*      t + 128 range is -128.0  +127.0
  55                      b*
  56
  57              16 bit Lab PCS:
  58
  59                      L*     0..100  into a 0..ff00 word.
  60                      a*     t + 128  range is  -128.0  +127.9961
  61                      b*
  62
  63
  64        We are always playing with 16 bits-data, so I will ignore the
  65        8-bits encoding scheme.
  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
  91
  92
  93 // On most modern computers, D > 4 M (i.e. a division takes more than 4
  94 // multiplications worth of time), so it is probably preferable to compute
  95 // a 24 bit result directly.
  96
  97 // #define ITERATE 1
  98
  99 static
 100 float CubeRoot(float x)
 101 {
 102        float fr, r;
 103        int ex, shx;
 104
 105        /* Argument reduction */
 106        fr = (float) frexp(x, &ex); /* separate into mantissa and exponent */
 107        shx = ex % 3;
 108
 109        if (shx > 0)
 110               shx -= 3; /* compute shx such that (ex - shx) is divisible by 3 */
 111
 112        ex = (ex - shx) / 3;        /* exponent of cube root */
 113        fr = (float) ldexp(fr, shx);
 114
 115        /* 0.125 <= fr < 1.0 */
 116
 117 #ifdef ITERATE
 118        /* Compute seed with a quadratic approximation */
 119
 120        fr = (-0.46946116F * fr + 1.072302F) * fr + 0.3812513F;/* 0.5<=fr<1 */
 121        r = ldexp(fr, ex);          /* 6 bits of precision */
 122
 123        /* Newton-Raphson iterations */
 124
 125        r = (float)(2.0/3.0) * r + (float)(1.0/3.0) * x / (r * r); /* 12 bits */
 126        r = (float)(2.0/3.0) * r + (float)(1.0/3.0) * x / (r * r); /* 24 bits */
 127 #else /* ITERATE */
 128
 129        /* Use quartic rational polynomial with error < 2^(-24) */
 130
 131        fr = (float) (((((45.2548339756803022511987494 * fr +
 132        192.2798368355061050458134625) * fr +
 133        119.1654824285581628956914143) * fr +
 134        13.43250139086239872172837314) * fr +
 135        0.1636161226585754240958355063)
 136        /
 137        ((((14.80884093219134573786480845 * fr +
 138        151.9714051044435648658557668) * fr +
 139        168.5254414101568283957668343) * fr +
 140        33.9905941350215598754191872) * fr +
 141        1.0));
 142        r = (float) ldexp(fr, ex); /* 24 bits of precision */
 143 #endif
 144        return r;
 145 }
 146
 147 static
 148 double f(double t)
 149 {
 150
 151       const double Limit = (24.0/116.0) * (24.0/116.0) * (24.0/116.0);
 152
 153        if (t <= Limit)
 154               return (841.0/108.0) * t + (16.0/116.0);
 155        else
 156               return CubeRoot((float) t);
 157 }
 158
 159
 160 static
 161 double f_1(double t)
 162 {
 163        const double Limit = (24.0/116.0);
 164
 165        if (t <= Limit)
 166        {
 167               double tmp;
 168
 169               tmp = (108.0/841.0) * (t - (16.0/116.0));
 170               if (tmp <= 0.0) return 0.0;
 171               else return tmp;
 172        }
 173
 174        return t * t * t;
 175 }
 176
 177
 178
 179 void LCMSEXPORT cmsXYZ2Lab(LPcmsCIEXYZ WhitePoint, LPcmsCIELab Lab, const cmsCIEXYZ* xyz)
 180 {
 181        double fx, fy, fz;
 182
 183        if (xyz -> X == 0 && xyz -> Y == 0 && xyz -> Z == 0)
 184        {
 185         Lab -> L = 0;
 186         Lab -> a = 0;
 187         Lab -> b = 0;
 188         return;
 189        }
 190
 191        if (WhitePoint == NULL)
 192             WhitePoint = cmsD50_XYZ();
 193
 194        fx = f(xyz->X / WhitePoint->X);
 195        fy = f(xyz->Y / WhitePoint->Y);
 196        fz = f(xyz->Z / WhitePoint->Z);
 197
 198        Lab->L = 116.0* fy - 16.;
 199
 200        Lab->a = 500.0*(fx - fy);
 201        Lab->b = 200.0*(fy - fz);
 202 }
 203
 204
 205
 206 void cmsXYZ2LabEncoded(WORD XYZ[3], WORD Lab[3])
 207 {
 208        Fixed32 X, Y, Z;
 209        double x, y, z, L, a, b;
 210        double fx, fy, fz;
 211        Fixed32 wL, wa, wb;
 212
 213        X = (Fixed32) XYZ[0] << 1;
 214        Y = (Fixed32) XYZ[1] << 1;
 215        Z = (Fixed32) XYZ[2] << 1;
 216
 217
 218        if (X==0 && Y==0 && Z==0) {
 219
 220                      Lab[0] = 0;
 221                      Lab[1] = Lab[2] =  0x8000;
 222                      return;
 223        }
 224
 225        // PCS is in D50
 226
 227
 228        x = FIXED_TO_DOUBLE(X) / D50X;
 229        y = FIXED_TO_DOUBLE(Y) / D50Y;
 230        z = FIXED_TO_DOUBLE(Z) / D50Z;
 231
 232
 233        fx = f(x);
 234        fy = f(y);
 235        fz = f(z);
 236
 237        L = 116.* fy - 16.;
 238
 239        a = 500.*(fx - fy);
 240        b = 200.*(fy - fz);
 241
 242        a += 128.;
 243        b += 128.;
 244
 245        wL = (int) (L * 652.800 + .5);
 246        wa = (int) (a * 256.0   + .5);
 247        wb = (int) (b * 256.0   + .5);
 248
 249
 250        Lab[0] = Clamp_L(wL);
 251        Lab[1] = Clamp_ab(wa);
 252        Lab[2] = Clamp_ab(wb);
 253
 254
 255 }
 256
 257
 258
 259
 260
 261
 262 void LCMSEXPORT cmsLab2XYZ(LPcmsCIEXYZ WhitePoint, LPcmsCIEXYZ xyz,  const cmsCIELab* Lab)
 263 {
 264         double x, y, z;
 265
 266         if (Lab -> L <= 0) {
 267                xyz -> X = 0;
 268                xyz -> Y = 0;
 269                xyz -> Z = 0;
 270                return;
 271         }
 272
 273
 274        if (WhitePoint == NULL)
 275             WhitePoint = cmsD50_XYZ();
 276
 277        y = (Lab-> L + 16.0) / 116.0;
 278        x = y + 0.002 * Lab -> a;
 279        z = y - 0.005 * Lab -> b;
 280
 281        xyz -> X = f_1(x) * WhitePoint -> X;
 282        xyz -> Y = f_1(y) * WhitePoint -> Y;
 283        xyz -> Z = f_1(z) * WhitePoint -> Z;
 284
 285 }
 286
 287
 288
 289 void cmsLab2XYZEncoded(WORD Lab[3], WORD XYZ[3])
 290 {
 291        double L, a, b;
 292        double X, Y, Z, x, y, z;
 293
 294
 295        L = ((double) Lab[0] * 100.0) / 65280.0;
 296        if (L==0.0) {
 297
 298        XYZ[0] = 0; XYZ[1] = 0; XYZ[2] = 0;
 299        return;
 300        }
 301
 302        a = ((double) Lab[1] / 256.0) - 128.0;
 303        b = ((double) Lab[2] / 256.0) - 128.0;
 304
 305        y = (L + 16.) / 116.0;
 306        x = y + 0.002 * a;
 307        z = y - 0.005 * b;
 308
 309        X = f_1(x) * D50X;
 310        Y = f_1(y) * D50Y;
 311        Z = f_1(z) * D50Z;
 312
 313        // Convert to 1.15 fixed format PCS
 314
 315
 316        XYZ[0] = _cmsClampWord((int) floor(X * 32768.0 + 0.5));
 317        XYZ[1] = _cmsClampWord((int) floor(Y * 32768.0 + 0.5));
 318        XYZ[2] = _cmsClampWord((int) floor(Z * 32768.0 + 0.5));
 319
 320
 321 }
 322
 323 static
 324 double L2float3(WORD v)
 325 {
 326        Fixed32 fix32;
 327
 328        fix32 = (Fixed32) v;
 329        return (double) fix32 / 652.800;
 330 }
 331
 332
 333 // the a/b part
 334
 335 static
 336 double ab2float3(WORD v)
 337 {
 338        Fixed32 fix32;
 339
 340        fix32 = (Fixed32) v;
 341        return ((double) fix32/256.0)-128.0;
 342 }
 343
 344 static
 345 WORD L2Fix3(double L)
 346 {
 347         return (WORD) (L *  652.800 + 0.5);
 348 }
 349
 350 static
 351 WORD ab2Fix3(double ab)
 352 {
 353         return (WORD) ((ab + 128.0) * 256.0 + 0.5);
 354 }
 355
 356
 357 // ICC 4.0 -- ICC has changed PCS Lab encoding.
 358
 359 static
 360 WORD L2Fix4(double L)
 361 {
 362      return (WORD) (L *  655.35 + 0.5);
 363 }
 364
 365 static
 366 WORD ab2Fix4(double ab)
 367 {
 368         return (WORD) ((ab + 128.0) * 257.0 + 0.5);
 369 }
 370
 371 static
 372 double L2float4(WORD v)
 373 {
 374        Fixed32 fix32;
 375
 376        fix32 = (Fixed32) v;
 377        return (double) fix32 / 655.35;
 378 }
 379
 380
 381 // the a/b part
 382
 383 static
 384 double ab2float4(WORD v)
 385 {
 386        Fixed32 fix32;
 387
 388        fix32 = (Fixed32) v;
 389        return ((double) fix32/257.0)-128.0;
 390 }
 391
 392
 393 void LCMSEXPORT cmsLabEncoded2Float(LPcmsCIELab Lab, const WORD wLab[3])
 394 {
 395         Lab->L = L2float3(wLab[0]);
 396         Lab->a = ab2float3(wLab[1]);
 397         Lab->b = ab2float3(wLab[2]);
 398 }
 399
 400
 401 void LCMSEXPORT cmsLabEncoded2Float4(LPcmsCIELab Lab, const WORD wLab[3])
 402 {
 403         Lab->L = L2float4(wLab[0]);
 404         Lab->a = ab2float4(wLab[1]);
 405         Lab->b = ab2float4(wLab[2]);
 406 }
 407
 408 static
 409 double Clamp_L_double(double L)
 410 {
 411     if (L < 0) L = 0;
 412     if (L > 100) L = 100;
 413
 414     return L;
 415 }
 416
 417
 418 static
 419 double Clamp_ab_double(double ab)
 420 {
 421     if (ab < -128) ab = -128.0;
 422     if (ab > +127.9961) ab = +127.9961;
 423
 424     return ab;
 425 }
 426
 427 void LCMSEXPORT cmsFloat2LabEncoded(WORD wLab[3], const cmsCIELab* fLab)
 428 {
 429     cmsCIELab Lab;
 430
 431
 432     Lab.L = Clamp_L_double(fLab ->L);
 433     Lab.a = Clamp_ab_double(fLab ->a);
 434     Lab.b = Clamp_ab_double(fLab ->b);
 435
 436     wLab[0] = L2Fix3(Lab.L);
 437     wLab[1] = ab2Fix3(Lab.a);
 438     wLab[2] = ab2Fix3(Lab.b);
 439 }
 440
 441
 442 void LCMSEXPORT cmsFloat2LabEncoded4(WORD wLab[3], const cmsCIELab* fLab)
 443 {
 444     cmsCIELab Lab;
 445
 446
 447     Lab.L = fLab ->L;
 448     Lab.a = fLab ->a;
 449     Lab.b = fLab ->b;
 450
 451
 452     if (Lab.L < 0) Lab.L = 0;
 453     if (Lab.L > 100.) Lab.L = 100.;
 454
 455     if (Lab.a < -128.) Lab.a = -128.;
 456     if (Lab.a > 127.) Lab.a = 127.;
 457     if (Lab.b < -128.) Lab.b = -128.;
 458     if (Lab.b > 127.) Lab.b = 127.;
 459
 460
 461     wLab[0] = L2Fix4(Lab.L);
 462     wLab[1] = ab2Fix4(Lab.a);
 463     wLab[2] = ab2Fix4(Lab.b);
 464 }
 465
 466
 467
 468
 469 void LCMSEXPORT cmsLab2LCh(LPcmsCIELCh LCh, const cmsCIELab* Lab)
 470 {
 471     double a, b;
 472
 473     LCh -> L = Clamp_L_double(Lab -> L);
 474
 475     a = Clamp_ab_double(Lab -> a);
 476     b = Clamp_ab_double(Lab -> b);
 477
 478     LCh -> C = pow(a * a + b * b, 0.5);
 479
 480     if (a == 0 && b == 0)
 481             LCh -> h   = 0;
 482     else
 483             LCh -> h = atan2(b, a);
 484
 485
 486     LCh -> h *= (180. / M_PI);
 487
 488
 489     while (LCh -> h >= 360.)         // Not necessary, but included as a check.
 490                 LCh -> h -= 360.;
 491
 492     while (LCh -> h < 0)
 493                 LCh -> h += 360.;
 494
 495 }
 496
 497
 498
 499
 500 void LCMSEXPORT cmsLCh2Lab(LPcmsCIELab Lab, const cmsCIELCh* LCh)
 501 {
 502
 503     double h = (LCh -> h * M_PI) / 180.0;
 504
 505     Lab -> L = Clamp_L_double(LCh -> L);
 506     Lab -> a = Clamp_ab_double(LCh -> C * cos(h));
 507     Lab -> b = Clamp_ab_double(LCh -> C * sin(h));
 508
 509 }
 510
 511
 512
 513
 514
 515 // In XYZ All 3 components are encoded using 1.15 fixed point
 516
 517 static
 518 WORD XYZ2Fix(double d)
 519 {
 520     return (WORD) floor(d * 32768.0 + 0.5);
 521 }
 522
 523
 524 void LCMSEXPORT cmsFloat2XYZEncoded(WORD XYZ[3], const cmsCIEXYZ* fXYZ)
 525 {
 526     cmsCIEXYZ xyz;
 527
 528     xyz.X = fXYZ -> X;
 529     xyz.Y = fXYZ -> Y;
 530     xyz.Z = fXYZ -> Z;
 531
 532
 533     // Clamp to encodeable values.
 534     // 1.99997 is reserved as out-of-gamut marker
 535
 536
 537     if (xyz.Y <= 0) {
 538
 539                 xyz.X = 0;
 540                 xyz.Y = 0;
 541                 xyz.Z = 0;
 542     }
 543
 544
 545     if (xyz.X > 1.99996)
 546            xyz.X = 1.99996;
 547
 548     if (xyz.X < 0)
 549            xyz.X = 0;
 550
 551     if (xyz.Y > 1.99996)
 552                 xyz.Y = 1.99996;
 553
 554     if (xyz.Y < 0)
 555            xyz.Y = 0;
 556
 557
 558     if (xyz.Z > 1.99996)
 559                 xyz.Z = 1.99996;
 560
 561     if (xyz.Z < 0)
 562            xyz.Z = 0;
 563
 564
 565
 566     XYZ[0] = XYZ2Fix(xyz.X);
 567     XYZ[1] = XYZ2Fix(xyz.Y);
 568     XYZ[2] = XYZ2Fix(xyz.Z);
 569
 570 }
 571
 572
 573 //  To convert from Fixed 1.15 point to double
 574
 575 static
 576 double XYZ2float(WORD v)
 577 {
 578        Fixed32 fix32;
 579
 580        // From 1.15 to 15.16
 581
 582        fix32 = v << 1;
 583
 584        // From fixed 15.16 to double
 585
 586        return FIXED_TO_DOUBLE(fix32);
 587 }
 588
 589
 590 void LCMSEXPORT cmsXYZEncoded2Float(LPcmsCIEXYZ fXYZ, const WORD XYZ[3])
 591 {
 592
 593     fXYZ -> X = XYZ2float(XYZ[0]);
 594     fXYZ -> Y = XYZ2float(XYZ[1]);
 595     fXYZ -> Z = XYZ2float(XYZ[2]);
 596
 597 }
 598
 599
 600
 601