modules/lcms/src/cmscam02.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
  24
  25 // CIECAM 02 appearance model. Many thanks to Jordi Vilar for the debugging.
  26
  27 #include "lcms.h"
  28
  29
  30 LCMSAPI LCMSHANDLE LCMSEXPORT cmsCIECAM02Init(LPcmsViewingConditions pVC);
  31 LCMSAPI void       LCMSEXPORT cmsCIECAM02Done(LCMSHANDLE hModel);
  32 LCMSAPI void       LCMSEXPORT cmsCIECAM02Forward(LCMSHANDLE hModel, LPcmsCIEXYZ pIn, LPcmsJCh pOut);
  33 LCMSAPI void       LCMSEXPORT cmsCIECAM02Reverse(LCMSHANDLE hModel, LPcmsJCh pIn,    LPcmsCIEXYZ pOut);
  34
  35
  36 // ---------- Implementation --------------------------------------------
  37
  38 typedef struct  {
  39
  40     double XYZ[3];
  41     double RGB[3];
  42     double RGBc[3];
  43     double RGBp[3];
  44     double RGBpa[3];
  45     double a, b, h, e, H, A, J, Q, s, t, C, M;
  46     double abC[2];
  47     double abs[2];
  48     double abM[2];
  49
  50 } CAM02COLOR, *LPCAM02COLOR;
  51
  52 typedef struct  {
  53
  54     CAM02COLOR adoptedWhite;
  55     double LA, Yb;
  56     double F, c, Nc;
  57     int surround;
  58     double n, Nbb, Ncb, z, FL, D;
  59
  60 } cmsCIECAM02, *LPcmsCIECAM02;
  61
  62
  63 static
  64 double compute_n(LPcmsCIECAM02 pMod)
  65 {
  66     return(pMod -> Yb / pMod -> adoptedWhite.XYZ[1]);
  67 }
  68
  69 static
  70 double compute_z(LPcmsCIECAM02 pMod)
  71 {
  72     return(1.48 + pow(pMod -> n, 0.5));
  73 }
  74
  75 static
  76 double computeNbb(LPcmsCIECAM02 pMod)
  77 {
  78     return(0.725 * pow((1.0 / pMod -> n), 0.2));
  79 }
  80
  81 static
  82 double computeFL(LPcmsCIECAM02 pMod)
  83 {
  84     double k, FL;
  85
  86     k = 1.0 / ((5.0 * pMod->LA) + 1.0);
  87     FL = 0.2 * pow(k, 4.0) * (5.0 * pMod->LA) + 0.1 *
  88         (pow((1.0 - pow(k, 4.0)), 2.0)) *
  89         (pow((5.0 * pMod->LA), (1.0 / 3.0)));
  90
  91     return FL;
  92 }
  93
  94 static
  95 double computeD(LPcmsCIECAM02 pMod)
  96 {
  97     double D;
  98
  99     D = pMod->F - (1.0/3.6)*(exp(((-pMod ->LA-42) / 92.0)));
 100
 101     return D;
 102 }
 103
 104
 105 static
 106 CAM02COLOR XYZtoCAT02(CAM02COLOR clr)
 107 {
 108     clr.RGB[0] = (clr.XYZ[0] *  0.7328) + (clr.XYZ[1] *  0.4296) + (clr.XYZ[2] * -0.1624);
 109     clr.RGB[1] = (clr.XYZ[0] * -0.7036) + (clr.XYZ[1] *  1.6975) + (clr.XYZ[2] *  0.0061);
 110     clr.RGB[2] = (clr.XYZ[0] *  0.0030) + (clr.XYZ[1] *  0.0136) + (clr.XYZ[2] *  0.9834);
 111
 112     return clr;
 113 }
 114
 115 static
 116 CAM02COLOR ChromaticAdaptation(CAM02COLOR clr, LPcmsCIECAM02 pMod)
 117 {
 118     int i;
 119     for (i = 0; i < 3; i++) {
 120         clr.RGBc[i] = ((pMod -> adoptedWhite.XYZ[1] *
 121             (pMod->D / pMod -> adoptedWhite.RGB[i])) +
 122             (1.0 - pMod->D)) * clr.RGB[i];
 123     }
 124
 125     return clr;
 126 }
 127
 128
 129 static
 130 CAM02COLOR CAT02toHPE (CAM02COLOR clr)
 131 {
 132
 133     double M[9];
 134
 135
 136     M[0] =(( 0.38971 *  1.096124) + (0.68898 * 0.454369) + (-0.07868 * -0.009628));
 137     M[1] =(( 0.38971 * -0.278869) + (0.68898 * 0.473533) + (-0.07868 * -0.005698));
 138     M[2] =(( 0.38971 *  0.182745) + (0.68898 * 0.072098) + (-0.07868 *  1.015326));
 139     M[3] =((-0.22981 *  1.096124) + (1.18340 * 0.454369) + ( 0.04641 * -0.009628));
 140     M[4] =((-0.22981 * -0.278869) + (1.18340 * 0.473533) + ( 0.04641 * -0.005698));
 141     M[5] =((-0.22981 *  0.182745) + (1.18340 * 0.072098) + ( 0.04641 *  1.015326));
 142     M[6] =(-0.009628);
 143     M[7] =(-0.005698);
 144     M[8] =( 1.015326);
 145
 146     clr.RGBp[0] = (clr.RGBc[0] * M[0]) +  (clr.RGBc[1] * M[1]) + (clr.RGBc[2] * M[2]);
 147     clr.RGBp[1] = (clr.RGBc[0] * M[3]) +  (clr.RGBc[1] * M[4]) + (clr.RGBc[2] * M[5]);
 148     clr.RGBp[2] = (clr.RGBc[0] * M[6]) +  (clr.RGBc[1] * M[7]) + (clr.RGBc[2] * M[8]);
 149
 150     return  clr;
 151 }
 152
 153 static
 154 CAM02COLOR NonlinearCompression(CAM02COLOR clr, LPcmsCIECAM02 pMod)
 155 {
 156     int i;
 157     double temp;
 158
 159     for (i = 0; i < 3; i++) {
 160         if (clr.RGBp[i] < 0) {
 161
 162             temp = pow((-1.0 * pMod->FL * clr.RGBp[i] / 100.0), 0.42);
 163             clr.RGBpa[i] = (-1.0 * 400.0 * temp) / (temp + 27.13) + 0.1;
 164         }
 165         else {
 166             temp = pow((pMod->FL * clr.RGBp[i] / 100.0), 0.42);
 167             clr.RGBpa[i] = (400.0 * temp) / (temp + 27.13) + 0.1;
 168         }
 169     }
 170
 171     clr.A = (((2.0 * clr.RGBpa[0]) + clr.RGBpa[1] +
 172         (clr.RGBpa[2] / 20.0)) - 0.305) * pMod->Nbb;
 173
 174     return clr;
 175 }
 176
 177 static
 178 CAM02COLOR ComputeCorrelates(CAM02COLOR clr, LPcmsCIECAM02 pMod)
 179 {
 180     double a, b, temp, e, t, r2d, d2r;
 181
 182     a = clr.RGBpa[0] - (12.0 * clr.RGBpa[1] / 11.0) + (clr.RGBpa[2] / 11.0);
 183     b = (clr.RGBpa[0] + clr.RGBpa[1] - (2.0 * clr.RGBpa[2])) / 9.0;
 184
 185     r2d = (180.0 / 3.141592654);
 186     if (a == 0) {
 187         if (b == 0)     clr.h = 0;
 188         else if (b > 0) clr.h = 90;
 189         else            clr.h = 270;
 190     }
 191     else if (a > 0) {
 192         temp = b / a;
 193         if (b > 0)       clr.h = (r2d * atan(temp));
 194         else if (b == 0) clr.h = 0;
 195         else             clr.h = (r2d * atan(temp)) + 360;
 196     }
 197     else {
 198         temp = b / a;
 199         clr.h = (r2d * atan(temp)) + 180;
 200     }
 201
 202     d2r = (3.141592654 / 180.0);
 203     e = ((12500.0 / 13.0) * pMod->Nc * pMod->Ncb) *
 204         (cos((clr.h * d2r + 2.0)) + 3.8);
 205
 206     if (clr.h < 20.14) {
 207         temp = ((clr.h + 122.47)/1.2) + ((20.14 - clr.h)/0.8);
 208         clr.H = 300 + (100*((clr.h + 122.47)/1.2)) / temp;
 209     }
 210     else if (clr.h < 90.0) {
 211         temp = ((clr.h - 20.14)/0.8) + ((90.00 - clr.h)/0.7);
 212         clr.H = (100*((clr.h - 20.14)/0.8)) / temp;
 213     }
 214     else if (clr.h < 164.25) {
 215         temp = ((clr.h - 90.00)/0.7) + ((164.25 - clr.h)/1.0);
 216         clr.H = 100 + ((100*((clr.h - 90.00)/0.7)) / temp);
 217     }
 218     else if (clr.h < 237.53) {
 219         temp = ((clr.h - 164.25)/1.0) + ((237.53 - clr.h)/1.2);
 220         clr.H = 200 + ((100*((clr.h - 164.25)/1.0)) / temp);
 221     }
 222     else {
 223         temp = ((clr.h - 237.53)/1.2) + ((360 - clr.h + 20.14)/0.8);
 224         clr.H = 300 + ((100*((clr.h - 237.53)/1.2)) / temp);
 225     }
 226
 227     clr.J = 100.0 * pow((clr.A / pMod->adoptedWhite.A),
 228         (pMod->c * pMod->z));
 229
 230     clr.Q = (4.0 / pMod->c) * pow((clr.J / 100.0), 0.5) *
 231         (pMod->adoptedWhite.A + 4.0) * pow(pMod->FL, 0.25);
 232
 233     t = (e * pow(((a * a) + (b * b)), 0.5)) /
 234         (clr.RGBpa[0] + clr.RGBpa[1] +
 235         ((21.0 / 20.0) * clr.RGBpa[2]));
 236
 237     clr.C = pow(t, 0.9) * pow((clr.J / 100.0), 0.5) *
 238         pow((1.64 - pow(0.29, pMod->n)), 0.73);
 239
 240     clr.M = clr.C * pow(pMod->FL, 0.25);
 241     clr.s = 100.0 * pow((clr.M / clr.Q), 0.5);
 242
 243     return clr;
 244 }
 245
 246
 247 static
 248 CAM02COLOR InverseCorrelates(CAM02COLOR clr, LPcmsCIECAM02 pMod)
 249 {
 250
 251     double t, e, p1, p2, p3, p4, p5, hr, d2r;
 252     d2r = 3.141592654 / 180.0;
 253
 254     t = pow( (clr.C / (pow((clr.J / 100.0), 0.5) *
 255         (pow((1.64 - pow(0.29, pMod->n)), 0.73)))),
 256         (1.0 / 0.9) );
 257     e = ((12500.0 / 13.0) * pMod->Nc * pMod->Ncb) *
 258         (cos((clr.h * d2r + 2.0)) + 3.8);
 259
 260     clr.A = pMod->adoptedWhite.A * pow(
 261            (clr.J / 100.0),
 262            (1.0 / (pMod->c * pMod->z)));
 263
 264     p1 = e / t;
 265     p2 = (clr.A / pMod->Nbb) + 0.305;
 266     p3 = 21.0 / 20.0;
 267
 268     hr = clr.h * d2r;
 269
 270     if (fabs(sin(hr)) >= fabs(cos(hr))) {
 271         p4 = p1 / sin(hr);
 272         clr.b = (p2 * (2.0 + p3) * (460.0 / 1403.0)) /
 273             (p4 + (2.0 + p3) * (220.0 / 1403.0) *
 274             (cos(hr) / sin(hr)) - (27.0 / 1403.0) +
 275             p3 * (6300.0 / 1403.0));
 276         clr.a = clr.b * (cos(hr) / sin(hr));
 277     }
 278     else {
 279         p5 = p1 / cos(hr);
 280         clr.a = (p2 * (2.0 + p3) * (460.0 / 1403.0)) /
 281             (p5 + (2.0 + p3) * (220.0 / 1403.0) -
 282             ((27.0 / 1403.0) - p3 * (6300.0 / 1403.0)) *
 283             (sin(hr) / cos(hr)));
 284         clr.b = clr.a * (sin(hr) / cos(hr));
 285     }
 286
 287     clr.RGBpa[0] = ((460.0 / 1403.0) * p2) +
 288               ((451.0 / 1403.0) * clr.a) +
 289               ((288.0 / 1403.0) * clr.b);
 290     clr.RGBpa[1] = ((460.0 / 1403.0) * p2) -
 291               ((891.0 / 1403.0) * clr.a) -
 292               ((261.0 / 1403.0) * clr.b);
 293     clr.RGBpa[2] = ((460.0 / 1403.0) * p2) -
 294               ((220.0 / 1403.0) * clr.a) -
 295               ((6300.0 / 1403.0) * clr.b);
 296
 297     return clr;
 298 }
 299
 300 static
 301 CAM02COLOR InverseNonlinearity(CAM02COLOR clr, LPcmsCIECAM02 pMod)
 302 {
 303     int i;
 304     double c1;
 305
 306     for (i = 0; i < 3; i++) {
 307         if ((clr.RGBpa[i] - 0.1) < 0) c1 = -1;
 308         else                               c1 = 1;
 309         clr.RGBp[i] = c1 * (100.0 / pMod->FL) *
 310             pow(((27.13 * fabs(clr.RGBpa[i] - 0.1)) /
 311             (400.0 - fabs(clr.RGBpa[i] - 0.1))),
 312             (1.0 / 0.42));
 313     }
 314
 315     return clr;
 316 }
 317
 318 static
 319 CAM02COLOR HPEtoCAT02(CAM02COLOR clr)
 320 {
 321     double M[9];
 322
 323     M[0] = (( 0.7328 *  1.910197) + (0.4296 * 0.370950));
 324     M[1] = (( 0.7328 * -1.112124) + (0.4296 * 0.629054));
 325     M[2] = (( 0.7328 *  0.201908) + (0.4296 * 0.000008) - 0.1624);
 326     M[3] = ((-0.7036 *  1.910197) + (1.6975 * 0.370950));
 327     M[4] = ((-0.7036 * -1.112124) + (1.6975 * 0.629054));
 328     M[5] = ((-0.7036 *  0.201908) + (1.6975 * 0.000008) + 0.0061);
 329     M[6] = (( 0.0030 *  1.910197) + (0.0136 * 0.370950));
 330     M[7] = (( 0.0030 * -1.112124) + (0.0136 * 0.629054));
 331     M[8] = (( 0.0030 *  0.201908) + (0.0136 * 0.000008) + 0.9834);;
 332
 333     clr.RGBc[0] = (clr.RGBp[0] * M[0]) + (clr.RGBp[1] * M[1]) + (clr.RGBp[2] * M[2]);
 334     clr.RGBc[1] = (clr.RGBp[0] * M[3]) + (clr.RGBp[1] * M[4]) + (clr.RGBp[2] * M[5]);
 335     clr.RGBc[2] = (clr.RGBp[0] * M[6]) + (clr.RGBp[1] * M[7]) + (clr.RGBp[2] * M[8]);
 336     return (clr);
 337 }
 338
 339
 340 static
 341 CAM02COLOR InverseChromaticAdaptation(CAM02COLOR clr,  LPcmsCIECAM02 pMod)
 342 {
 343     int i;
 344     for (i = 0; i < 3; i++) {
 345         clr.RGB[i] = clr.RGBc[i] /
 346             ((pMod->adoptedWhite.XYZ[1] * pMod->D / pMod->adoptedWhite.RGB[i]) + 1.0 - pMod->D);
 347     }
 348     return(clr);
 349 }
 350
 351
 352 static
 353 CAM02COLOR CAT02toXYZ(CAM02COLOR clr)
 354 {
 355     clr.XYZ[0] = (clr.RGB[0] *  1.096124) + (clr.RGB[1] * -0.278869) + (clr.RGB[2] *  0.182745);
 356     clr.XYZ[1] = (clr.RGB[0] *  0.454369) + (clr.RGB[1] *  0.473533) + (clr.RGB[2] *  0.072098);
 357     clr.XYZ[2] = (clr.RGB[0] * -0.009628) + (clr.RGB[1] * -0.005698) + (clr.RGB[2] *  1.015326);
 358
 359     return(clr);
 360 }
 361
 362
 363
 364
 365 LCMSHANDLE LCMSEXPORT cmsCIECAM02Init(LPcmsViewingConditions pVC)
 366 {
 367     LPcmsCIECAM02 lpMod;
 368
 369
 370    if((lpMod = (LPcmsCIECAM02) _cmsMalloc(sizeof(cmsCIECAM02))) == NULL) {
 371         return (LCMSHANDLE) NULL;
 372     }
 373
 374
 375     ZeroMemory(lpMod, sizeof(cmsCIECAM02));
 376
 377     lpMod ->adoptedWhite.XYZ[0] = pVC ->whitePoint.X;
 378     lpMod ->adoptedWhite.XYZ[1] = pVC ->whitePoint.Y;
 379     lpMod ->adoptedWhite.XYZ[2] = pVC ->whitePoint.Z;
 380
 381     lpMod -> LA       = pVC ->La;
 382     lpMod -> Yb       = pVC ->Yb;
 383     lpMod -> D        = pVC ->D_value;
 384     lpMod -> surround = pVC ->surround;
 385
 386     switch (lpMod -> surround) {
 387
 388     case AVG_SURROUND_4:
 389         lpMod->F = 1.0;     // Not included in CAM02
 390         lpMod->c = 0.69;
 391         lpMod->Nc = 1.0;
 392         break;
 393
 394     case CUTSHEET_SURROUND:
 395         lpMod->F = 0.8;
 396         lpMod->c = 0.41;
 397         lpMod->Nc = 0.8;
 398         break;
 399
 400     case DARK_SURROUND:
 401         lpMod -> F  = 0.8;
 402         lpMod -> c  = 0.525;
 403         lpMod -> Nc = 0.8;
 404         break;
 405
 406
 407     case DIM_SURROUND:
 408         lpMod -> F  = 0.9;
 409         lpMod -> c  = 0.59;
 410         lpMod -> Nc = 0.95;
 411         break;
 412
 413     default:
 414         // Average surround
 415         lpMod -> F  = 1.0;
 416         lpMod -> c  = 0.69;
 417         lpMod -> Nc = 1.0;
 418     }
 419
 420     lpMod -> n   = compute_n(lpMod);
 421     lpMod -> z   = compute_z(lpMod);
 422     lpMod -> Nbb = computeNbb(lpMod);
 423     lpMod -> FL  = computeFL(lpMod);
 424
 425     if (lpMod -> D == D_CALCULATE ||
 426         lpMod -> D == D_CALCULATE_DISCOUNT) {
 427
 428     lpMod -> D   = computeD(lpMod);
 429     }
 430
 431     lpMod -> Ncb = lpMod -> Nbb;
 432
 433     lpMod -> adoptedWhite = XYZtoCAT02(lpMod -> adoptedWhite);
 434     lpMod -> adoptedWhite = ChromaticAdaptation(lpMod -> adoptedWhite, lpMod);
 435     lpMod -> adoptedWhite = CAT02toHPE(lpMod -> adoptedWhite);
 436     lpMod -> adoptedWhite = NonlinearCompression(lpMod -> adoptedWhite, lpMod);
 437
 438     return (LCMSHANDLE) lpMod;
 439
 440 }
 441
 442 void LCMSEXPORT cmsCIECAM02Done(LCMSHANDLE hModel)
 443 {
 444     LPcmsCIECAM02 lpMod = (LPcmsCIECAM02) (LPSTR) hModel;
 445     if (lpMod) _cmsFree(lpMod);
 446 }
 447
 448
 449 void LCMSEXPORT cmsCIECAM02Forward(LCMSHANDLE hModel, LPcmsCIEXYZ pIn, LPcmsJCh pOut)
 450 {
 451     CAM02COLOR clr;
 452     LPcmsCIECAM02 lpMod = (LPcmsCIECAM02) (LPSTR) hModel;
 453
 454     clr.XYZ[0] = pIn ->X;
 455     clr.XYZ[1] = pIn ->Y;
 456     clr.XYZ[2] = pIn ->Z;
 457
 458     clr = XYZtoCAT02(clr);
 459     clr = ChromaticAdaptation(clr, lpMod);
 460     clr = CAT02toHPE(clr);
 461     clr = NonlinearCompression(clr, lpMod);
 462     clr = ComputeCorrelates(clr, lpMod);
 463
 464     pOut ->J = clr.J;
 465     pOut ->C = clr.C;
 466     pOut ->h = clr.h;
 467 }
 468
 469 void LCMSEXPORT cmsCIECAM02Reverse(LCMSHANDLE hModel, LPcmsJCh pIn, LPcmsCIEXYZ pOut)
 470 {
 471     CAM02COLOR clr;
 472     LPcmsCIECAM02 lpMod = (LPcmsCIECAM02) (LPSTR) hModel;
 473
 474
 475     clr.J = pIn -> J;
 476     clr.C = pIn -> C;
 477     clr.h = pIn -> h;
 478
 479     clr = InverseCorrelates(clr, lpMod);
 480     clr = InverseNonlinearity(clr, lpMod);
 481     clr = HPEtoCAT02(clr);
 482     clr = InverseChromaticAdaptation(clr, lpMod);
 483     clr = CAT02toXYZ(clr);
 484
 485     pOut ->X = clr.XYZ[0];
 486     pOut ->Y = clr.XYZ[1];
 487     pOut ->Z = clr.XYZ[2];
 488
 489 }
 490