1 /*
   2  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   3  *
   4  * This code is free software; you can redistribute it and/or modify it
   5  * under the terms of the GNU General Public License version 2 only, as
   6  * published by the Free Software Foundation.  Oracle designates this
   7  * particular file as subject to the "Classpath" exception as provided
   8  * by Oracle in the LICENSE file that accompanied this code.
   9  *
  10  * This code is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13  * version 2 for more details (a copy is included in the LICENSE file that
  14  * accompanied this code).
  15  *
  16  * You should have received a copy of the GNU General Public License version
  17  * 2 along with this work; if not, write to the Free Software Foundation,
  18  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
  19  *
  20  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
  21  * or visit www.oracle.com if you need additional information or have any
  22  * questions.
  23  */
  24 
  25 // This file is available under and governed by the GNU General Public
  26 // License version 2 only, as published by the Free Software Foundation.
  27 // However, the following notice accompanied the original version of this
  28 // file:
  29 //
  30 //---------------------------------------------------------------------------------
  31 //
  32 //  Little Color Management System
  33 //  Copyright (c) 1998-2020 Marti Maria Saguer
  34 //
  35 // Permission is hereby granted, free of charge, to any person obtaining
  36 // a copy of this software and associated documentation files (the "Software"),
  37 // to deal in the Software without restriction, including without limitation
  38 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  39 // and/or sell copies of the Software, and to permit persons to whom the Software
  40 // is furnished to do so, subject to the following conditions:
  41 //
  42 // The above copyright notice and this permission notice shall be included in
  43 // all copies or substantial portions of the Software.
  44 //
  45 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  46 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
  47 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  48 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  49 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  50 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  51 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  52 //
  53 //---------------------------------------------------------------------------------
  54 //
  55 
  56 #include "lcms2_internal.h"
  57 
  58 // CIECAM 02 appearance model. Many thanks to Jordi Vilar for the debugging.
  59 
  60 // ---------- Implementation --------------------------------------------
  61 
  62 typedef struct  {
  63 
  64     cmsFloat64Number XYZ[3];
  65     cmsFloat64Number RGB[3];
  66     cmsFloat64Number RGBc[3];
  67     cmsFloat64Number RGBp[3];
  68     cmsFloat64Number RGBpa[3];
  69     cmsFloat64Number a, b, h, e, H, A, J, Q, s, t, C, M;
  70     cmsFloat64Number abC[2];
  71     cmsFloat64Number abs[2];
  72     cmsFloat64Number abM[2];
  73 
  74 } CAM02COLOR;
  75 
  76 typedef struct  {
  77 
  78     CAM02COLOR adoptedWhite;
  79     cmsFloat64Number LA, Yb;
  80     cmsFloat64Number F, c, Nc;
  81     cmsUInt32Number surround;
  82     cmsFloat64Number n, Nbb, Ncb, z, FL, D;
  83 
  84     cmsContext ContextID;
  85 
  86 } cmsCIECAM02;
  87 
  88 
  89 static
  90 cmsFloat64Number compute_n(cmsCIECAM02* pMod)
  91 {
  92     return (pMod -> Yb / pMod -> adoptedWhite.XYZ[1]);
  93 }
  94 
  95 static
  96 cmsFloat64Number compute_z(cmsCIECAM02* pMod)
  97 {
  98     return (1.48 + pow(pMod -> n, 0.5));
  99 }
 100 
 101 static
 102 cmsFloat64Number computeNbb(cmsCIECAM02* pMod)
 103 {
 104     return (0.725 * pow((1.0 / pMod -> n), 0.2));
 105 }
 106 
 107 static
 108 cmsFloat64Number computeFL(cmsCIECAM02* pMod)
 109 {
 110     cmsFloat64Number k, FL;
 111 
 112     k = 1.0 / ((5.0 * pMod->LA) + 1.0);
 113     FL = 0.2 * pow(k, 4.0) * (5.0 * pMod->LA) + 0.1 *
 114         (pow((1.0 - pow(k, 4.0)), 2.0)) *
 115         (pow((5.0 * pMod->LA), (1.0 / 3.0)));
 116 
 117     return FL;
 118 }
 119 
 120 static
 121 cmsFloat64Number computeD(cmsCIECAM02* pMod)
 122 {
 123     cmsFloat64Number D;
 124 
 125     D = pMod->F - (1.0/3.6)*(exp(((-pMod ->LA-42) / 92.0)));
 126 
 127     return D;
 128 }
 129 
 130 
 131 static
 132 CAM02COLOR XYZtoCAT02(CAM02COLOR clr)
 133 {
 134     clr.RGB[0] = (clr.XYZ[0] *  0.7328) + (clr.XYZ[1] *  0.4296) + (clr.XYZ[2] * -0.1624);
 135     clr.RGB[1] = (clr.XYZ[0] * -0.7036) + (clr.XYZ[1] *  1.6975) + (clr.XYZ[2] *  0.0061);
 136     clr.RGB[2] = (clr.XYZ[0] *  0.0030) + (clr.XYZ[1] *  0.0136) + (clr.XYZ[2] *  0.9834);
 137 
 138     return clr;
 139 }
 140 
 141 static
 142 CAM02COLOR ChromaticAdaptation(CAM02COLOR clr, cmsCIECAM02* pMod)
 143 {
 144     cmsUInt32Number i;
 145 
 146     for (i = 0; i < 3; i++) {
 147         clr.RGBc[i] = ((pMod -> adoptedWhite.XYZ[1] *
 148             (pMod->D / pMod -> adoptedWhite.RGB[i])) +
 149             (1.0 - pMod->D)) * clr.RGB[i];
 150     }
 151 
 152     return clr;
 153 }
 154 
 155 
 156 static
 157 CAM02COLOR CAT02toHPE(CAM02COLOR clr)
 158 {
 159     cmsFloat64Number M[9];
 160 
 161     M[0] =(( 0.38971 *  1.096124) + (0.68898 * 0.454369) + (-0.07868 * -0.009628));
 162     M[1] =(( 0.38971 * -0.278869) + (0.68898 * 0.473533) + (-0.07868 * -0.005698));
 163     M[2] =(( 0.38971 *  0.182745) + (0.68898 * 0.072098) + (-0.07868 *  1.015326));
 164     M[3] =((-0.22981 *  1.096124) + (1.18340 * 0.454369) + ( 0.04641 * -0.009628));
 165     M[4] =((-0.22981 * -0.278869) + (1.18340 * 0.473533) + ( 0.04641 * -0.005698));
 166     M[5] =((-0.22981 *  0.182745) + (1.18340 * 0.072098) + ( 0.04641 *  1.015326));
 167     M[6] =(-0.009628);
 168     M[7] =(-0.005698);
 169     M[8] =( 1.015326);
 170 
 171     clr.RGBp[0] = (clr.RGBc[0] * M[0]) +  (clr.RGBc[1] * M[1]) + (clr.RGBc[2] * M[2]);
 172     clr.RGBp[1] = (clr.RGBc[0] * M[3]) +  (clr.RGBc[1] * M[4]) + (clr.RGBc[2] * M[5]);
 173     clr.RGBp[2] = (clr.RGBc[0] * M[6]) +  (clr.RGBc[1] * M[7]) + (clr.RGBc[2] * M[8]);
 174 
 175     return  clr;
 176 }
 177 
 178 static
 179 CAM02COLOR NonlinearCompression(CAM02COLOR clr, cmsCIECAM02* pMod)
 180 {
 181     cmsUInt32Number i;
 182     cmsFloat64Number temp;
 183 
 184     for (i = 0; i < 3; i++) {
 185         if (clr.RGBp[i] < 0) {
 186 
 187             temp = pow((-1.0 * pMod->FL * clr.RGBp[i] / 100.0), 0.42);
 188             clr.RGBpa[i] = (-1.0 * 400.0 * temp) / (temp + 27.13) + 0.1;
 189         }
 190         else {
 191             temp = pow((pMod->FL * clr.RGBp[i] / 100.0), 0.42);
 192             clr.RGBpa[i] = (400.0 * temp) / (temp + 27.13) + 0.1;
 193         }
 194     }
 195 
 196     clr.A = (((2.0 * clr.RGBpa[0]) + clr.RGBpa[1] +
 197         (clr.RGBpa[2] / 20.0)) - 0.305) * pMod->Nbb;
 198 
 199     return clr;
 200 }
 201 
 202 static
 203 CAM02COLOR ComputeCorrelates(CAM02COLOR clr, cmsCIECAM02* pMod)
 204 {
 205     cmsFloat64Number a, b, temp, e, t, r2d, d2r;
 206 
 207     a = clr.RGBpa[0] - (12.0 * clr.RGBpa[1] / 11.0) + (clr.RGBpa[2] / 11.0);
 208     b = (clr.RGBpa[0] + clr.RGBpa[1] - (2.0 * clr.RGBpa[2])) / 9.0;
 209 
 210     r2d = (180.0 / 3.141592654);
 211     if (a == 0) {
 212         if (b == 0)     clr.h = 0;
 213         else if (b > 0) clr.h = 90;
 214         else            clr.h = 270;
 215     }
 216     else if (a > 0) {
 217         temp = b / a;
 218         if (b > 0)       clr.h = (r2d * atan(temp));
 219         else if (b == 0) clr.h = 0;
 220         else             clr.h = (r2d * atan(temp)) + 360;
 221     }
 222     else {
 223         temp = b / a;
 224         clr.h = (r2d * atan(temp)) + 180;
 225     }
 226 
 227     d2r = (3.141592654 / 180.0);
 228     e = ((12500.0 / 13.0) * pMod->Nc * pMod->Ncb) *
 229         (cos((clr.h * d2r + 2.0)) + 3.8);
 230 
 231     if (clr.h < 20.14) {
 232         temp = ((clr.h + 122.47)/1.2) + ((20.14 - clr.h)/0.8);
 233         clr.H = 300 + (100*((clr.h + 122.47)/1.2)) / temp;
 234     }
 235     else if (clr.h < 90.0) {
 236         temp = ((clr.h - 20.14)/0.8) + ((90.00 - clr.h)/0.7);
 237         clr.H = (100*((clr.h - 20.14)/0.8)) / temp;
 238     }
 239     else if (clr.h < 164.25) {
 240         temp = ((clr.h - 90.00)/0.7) + ((164.25 - clr.h)/1.0);
 241         clr.H = 100 + ((100*((clr.h - 90.00)/0.7)) / temp);
 242     }
 243     else if (clr.h < 237.53) {
 244         temp = ((clr.h - 164.25)/1.0) + ((237.53 - clr.h)/1.2);
 245         clr.H = 200 + ((100*((clr.h - 164.25)/1.0)) / temp);
 246     }
 247     else {
 248         temp = ((clr.h - 237.53)/1.2) + ((360 - clr.h + 20.14)/0.8);
 249         clr.H = 300 + ((100*((clr.h - 237.53)/1.2)) / temp);
 250     }
 251 
 252     clr.J = 100.0 * pow((clr.A / pMod->adoptedWhite.A),
 253         (pMod->c * pMod->z));
 254 
 255     clr.Q = (4.0 / pMod->c) * pow((clr.J / 100.0), 0.5) *
 256         (pMod->adoptedWhite.A + 4.0) * pow(pMod->FL, 0.25);
 257 
 258     t = (e * pow(((a * a) + (b * b)), 0.5)) /
 259         (clr.RGBpa[0] + clr.RGBpa[1] +
 260         ((21.0 / 20.0) * clr.RGBpa[2]));
 261 
 262     clr.C = pow(t, 0.9) * pow((clr.J / 100.0), 0.5) *
 263         pow((1.64 - pow(0.29, pMod->n)), 0.73);
 264 
 265     clr.M = clr.C * pow(pMod->FL, 0.25);
 266     clr.s = 100.0 * pow((clr.M / clr.Q), 0.5);
 267 
 268     return clr;
 269 }
 270 
 271 
 272 static
 273 CAM02COLOR InverseCorrelates(CAM02COLOR clr, cmsCIECAM02* pMod)
 274 {
 275 
 276     cmsFloat64Number t, e, p1, p2, p3, p4, p5, hr, d2r;
 277     d2r = 3.141592654 / 180.0;
 278 
 279     t = pow( (clr.C / (pow((clr.J / 100.0), 0.5) *
 280         (pow((1.64 - pow(0.29, pMod->n)), 0.73)))),
 281         (1.0 / 0.9) );
 282     e = ((12500.0 / 13.0) * pMod->Nc * pMod->Ncb) *
 283         (cos((clr.h * d2r + 2.0)) + 3.8);
 284 
 285     clr.A = pMod->adoptedWhite.A * pow(
 286            (clr.J / 100.0),
 287            (1.0 / (pMod->c * pMod->z)));
 288 
 289     p1 = e / t;
 290     p2 = (clr.A / pMod->Nbb) + 0.305;
 291     p3 = 21.0 / 20.0;
 292 
 293     hr = clr.h * d2r;
 294 
 295     if (fabs(sin(hr)) >= fabs(cos(hr))) {
 296         p4 = p1 / sin(hr);
 297         clr.b = (p2 * (2.0 + p3) * (460.0 / 1403.0)) /
 298             (p4 + (2.0 + p3) * (220.0 / 1403.0) *
 299             (cos(hr) / sin(hr)) - (27.0 / 1403.0) +
 300             p3 * (6300.0 / 1403.0));
 301         clr.a = clr.b * (cos(hr) / sin(hr));
 302     }
 303     else {
 304         p5 = p1 / cos(hr);
 305         clr.a = (p2 * (2.0 + p3) * (460.0 / 1403.0)) /
 306             (p5 + (2.0 + p3) * (220.0 / 1403.0) -
 307             ((27.0 / 1403.0) - p3 * (6300.0 / 1403.0)) *
 308             (sin(hr) / cos(hr)));
 309         clr.b = clr.a * (sin(hr) / cos(hr));
 310     }
 311 
 312     clr.RGBpa[0] = ((460.0 / 1403.0) * p2) +
 313               ((451.0 / 1403.0) * clr.a) +
 314               ((288.0 / 1403.0) * clr.b);
 315     clr.RGBpa[1] = ((460.0 / 1403.0) * p2) -
 316               ((891.0 / 1403.0) * clr.a) -
 317               ((261.0 / 1403.0) * clr.b);
 318     clr.RGBpa[2] = ((460.0 / 1403.0) * p2) -
 319               ((220.0 / 1403.0) * clr.a) -
 320               ((6300.0 / 1403.0) * clr.b);
 321 
 322     return clr;
 323 }
 324 
 325 static
 326 CAM02COLOR InverseNonlinearity(CAM02COLOR clr, cmsCIECAM02* pMod)
 327 {
 328     cmsUInt32Number i;
 329     cmsFloat64Number c1;
 330 
 331     for (i = 0; i < 3; i++) {
 332         if ((clr.RGBpa[i] - 0.1) < 0) c1 = -1;
 333         else                               c1 = 1;
 334         clr.RGBp[i] = c1 * (100.0 / pMod->FL) *
 335             pow(((27.13 * fabs(clr.RGBpa[i] - 0.1)) /
 336             (400.0 - fabs(clr.RGBpa[i] - 0.1))),
 337             (1.0 / 0.42));
 338     }
 339 
 340     return clr;
 341 }
 342 
 343 static
 344 CAM02COLOR HPEtoCAT02(CAM02COLOR clr)
 345 {
 346     cmsFloat64Number M[9];
 347 
 348     M[0] = (( 0.7328 *  1.910197) + (0.4296 * 0.370950));
 349     M[1] = (( 0.7328 * -1.112124) + (0.4296 * 0.629054));
 350     M[2] = (( 0.7328 *  0.201908) + (0.4296 * 0.000008) - 0.1624);
 351     M[3] = ((-0.7036 *  1.910197) + (1.6975 * 0.370950));
 352     M[4] = ((-0.7036 * -1.112124) + (1.6975 * 0.629054));
 353     M[5] = ((-0.7036 *  0.201908) + (1.6975 * 0.000008) + 0.0061);
 354     M[6] = (( 0.0030 *  1.910197) + (0.0136 * 0.370950));
 355     M[7] = (( 0.0030 * -1.112124) + (0.0136 * 0.629054));
 356     M[8] = (( 0.0030 *  0.201908) + (0.0136 * 0.000008) + 0.9834);;
 357 
 358     clr.RGBc[0] = (clr.RGBp[0] * M[0]) + (clr.RGBp[1] * M[1]) + (clr.RGBp[2] * M[2]);
 359     clr.RGBc[1] = (clr.RGBp[0] * M[3]) + (clr.RGBp[1] * M[4]) + (clr.RGBp[2] * M[5]);
 360     clr.RGBc[2] = (clr.RGBp[0] * M[6]) + (clr.RGBp[1] * M[7]) + (clr.RGBp[2] * M[8]);
 361     return clr;
 362 }
 363 
 364 
 365 static
 366 CAM02COLOR InverseChromaticAdaptation(CAM02COLOR clr,  cmsCIECAM02* pMod)
 367 {
 368     cmsUInt32Number i;
 369     for (i = 0; i < 3; i++) {
 370         clr.RGB[i] = clr.RGBc[i] /
 371             ((pMod->adoptedWhite.XYZ[1] * pMod->D / pMod->adoptedWhite.RGB[i]) + 1.0 - pMod->D);
 372     }
 373     return clr;
 374 }
 375 
 376 
 377 static
 378 CAM02COLOR CAT02toXYZ(CAM02COLOR clr)
 379 {
 380     clr.XYZ[0] = (clr.RGB[0] *  1.096124) + (clr.RGB[1] * -0.278869) + (clr.RGB[2] *  0.182745);
 381     clr.XYZ[1] = (clr.RGB[0] *  0.454369) + (clr.RGB[1] *  0.473533) + (clr.RGB[2] *  0.072098);
 382     clr.XYZ[2] = (clr.RGB[0] * -0.009628) + (clr.RGB[1] * -0.005698) + (clr.RGB[2] *  1.015326);
 383 
 384     return clr;
 385 }
 386 
 387 
 388 cmsHANDLE  CMSEXPORT cmsCIECAM02Init(cmsContext ContextID, const cmsViewingConditions* pVC)
 389 {
 390     cmsCIECAM02* lpMod;
 391 
 392     _cmsAssert(pVC != NULL);
 393 
 394     if((lpMod = (cmsCIECAM02*) _cmsMallocZero(ContextID, sizeof(cmsCIECAM02))) == NULL) {
 395         return NULL;
 396     }
 397 
 398     lpMod ->ContextID = ContextID;
 399 
 400     lpMod ->adoptedWhite.XYZ[0] = pVC ->whitePoint.X;
 401     lpMod ->adoptedWhite.XYZ[1] = pVC ->whitePoint.Y;
 402     lpMod ->adoptedWhite.XYZ[2] = pVC ->whitePoint.Z;
 403 
 404     lpMod -> LA       = pVC ->La;
 405     lpMod -> Yb       = pVC ->Yb;
 406     lpMod -> D        = pVC ->D_value;
 407     lpMod -> surround = pVC ->surround;
 408 
 409     switch (lpMod -> surround) {
 410 
 411 
 412     case CUTSHEET_SURROUND:
 413         lpMod->F = 0.8;
 414         lpMod->c = 0.41;
 415         lpMod->Nc = 0.8;
 416         break;
 417 
 418     case DARK_SURROUND:
 419         lpMod -> F  = 0.8;
 420         lpMod -> c  = 0.525;
 421         lpMod -> Nc = 0.8;
 422         break;
 423 
 424     case DIM_SURROUND:
 425         lpMod -> F  = 0.9;
 426         lpMod -> c  = 0.59;
 427         lpMod -> Nc = 0.95;
 428         break;
 429 
 430     default:
 431         // Average surround
 432         lpMod -> F  = 1.0;
 433         lpMod -> c  = 0.69;
 434         lpMod -> Nc = 1.0;
 435     }
 436 
 437     lpMod -> n   = compute_n(lpMod);
 438     lpMod -> z   = compute_z(lpMod);
 439     lpMod -> Nbb = computeNbb(lpMod);
 440     lpMod -> FL  = computeFL(lpMod);
 441 
 442     if (lpMod -> D == D_CALCULATE) {
 443         lpMod -> D   = computeD(lpMod);
 444     }
 445 
 446     lpMod -> Ncb = lpMod -> Nbb;
 447 
 448     lpMod -> adoptedWhite = XYZtoCAT02(lpMod -> adoptedWhite);
 449     lpMod -> adoptedWhite = ChromaticAdaptation(lpMod -> adoptedWhite, lpMod);
 450     lpMod -> adoptedWhite = CAT02toHPE(lpMod -> adoptedWhite);
 451     lpMod -> adoptedWhite = NonlinearCompression(lpMod -> adoptedWhite, lpMod);
 452 
 453     return (cmsHANDLE) lpMod;
 454 
 455 }
 456 
 457 void CMSEXPORT cmsCIECAM02Done(cmsHANDLE hModel)
 458 {
 459     cmsCIECAM02* lpMod = (cmsCIECAM02*) hModel;
 460 
 461     if (lpMod) _cmsFree(lpMod ->ContextID, lpMod);
 462 }
 463 
 464 
 465 void CMSEXPORT cmsCIECAM02Forward(cmsHANDLE hModel, const cmsCIEXYZ* pIn, cmsJCh* pOut)
 466 {
 467     CAM02COLOR clr;
 468     cmsCIECAM02* lpMod = (cmsCIECAM02*) hModel;
 469 
 470     _cmsAssert(lpMod != NULL);
 471     _cmsAssert(pIn != NULL);
 472     _cmsAssert(pOut != NULL);
 473 
 474     memset(&clr, 0, sizeof(clr));
 475 
 476     clr.XYZ[0] = pIn ->X;
 477     clr.XYZ[1] = pIn ->Y;
 478     clr.XYZ[2] = pIn ->Z;
 479 
 480     clr = XYZtoCAT02(clr);
 481     clr = ChromaticAdaptation(clr, lpMod);
 482     clr = CAT02toHPE(clr);
 483     clr = NonlinearCompression(clr, lpMod);
 484     clr = ComputeCorrelates(clr, lpMod);
 485 
 486     pOut ->J = clr.J;
 487     pOut ->C = clr.C;
 488     pOut ->h = clr.h;
 489 }
 490 
 491 void CMSEXPORT cmsCIECAM02Reverse(cmsHANDLE hModel, const cmsJCh* pIn, cmsCIEXYZ* pOut)
 492 {
 493     CAM02COLOR clr;
 494     cmsCIECAM02* lpMod = (cmsCIECAM02*) hModel;
 495 
 496     _cmsAssert(lpMod != NULL);
 497     _cmsAssert(pIn != NULL);
 498     _cmsAssert(pOut != NULL);
 499 
 500     memset(&clr, 0, sizeof(clr));
 501 
 502     clr.J = pIn -> J;
 503     clr.C = pIn -> C;
 504     clr.h = pIn -> h;
 505 
 506     clr = InverseCorrelates(clr, lpMod);
 507     clr = InverseNonlinearity(clr, lpMod);
 508     clr = HPEtoCAT02(clr);
 509     clr = InverseChromaticAdaptation(clr, lpMod);
 510     clr = CAT02toXYZ(clr);
 511 
 512     pOut ->X = clr.XYZ[0];
 513     pOut ->Y = clr.XYZ[1];
 514     pOut ->Z = clr.XYZ[2];
 515 }