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
compute_n(LPcmsCIECAM02 pMod)64 double compute_n(LPcmsCIECAM02 pMod)
65 {
66 return(pMod -> Yb / pMod -> adoptedWhite.XYZ[1]);
67 }
68
69 static
compute_z(LPcmsCIECAM02 pMod)70 double compute_z(LPcmsCIECAM02 pMod)
71 {
72 return(1.48 + pow(pMod -> n, 0.5));
73 }
74
75 static
computeNbb(LPcmsCIECAM02 pMod)76 double computeNbb(LPcmsCIECAM02 pMod)
77 {
78 return(0.725 * pow((1.0 / pMod -> n), 0.2));
79 }
80
81 static
computeFL(LPcmsCIECAM02 pMod)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
computeD(LPcmsCIECAM02 pMod)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
XYZtoCAT02(CAM02COLOR clr)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
ChromaticAdaptation(CAM02COLOR clr,LPcmsCIECAM02 pMod)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
CAT02toHPE(CAM02COLOR clr)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
NonlinearCompression(CAM02COLOR clr,LPcmsCIECAM02 pMod)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
ComputeCorrelates(CAM02COLOR clr,LPcmsCIECAM02 pMod)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
InverseCorrelates(CAM02COLOR clr,LPcmsCIECAM02 pMod)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
InverseNonlinearity(CAM02COLOR clr,LPcmsCIECAM02 pMod)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
HPEtoCAT02(CAM02COLOR clr)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
InverseChromaticAdaptation(CAM02COLOR clr,LPcmsCIECAM02 pMod)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
CAT02toXYZ(CAM02COLOR clr)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
cmsCIECAM02Init(LPcmsViewingConditions pVC)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
cmsCIECAM02Done(LCMSHANDLE hModel)442 void LCMSEXPORT cmsCIECAM02Done(LCMSHANDLE hModel)
443 {
444 LPcmsCIECAM02 lpMod = (LPcmsCIECAM02) (LPSTR) hModel;
445 if (lpMod) _cmsFree(lpMod);
446 }
447
448
cmsCIECAM02Forward(LCMSHANDLE hModel,LPcmsCIEXYZ pIn,LPcmsJCh pOut)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
cmsCIECAM02Reverse(LCMSHANDLE hModel,LPcmsJCh pIn,LPcmsCIEXYZ pOut)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
491