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