1
2 /*
3 * International Color Consortium color transform expanded support
4 *
5 * Author: Graeme W. Gill
6 * Date: 2/7/00
7 * Version: 1.00
8 *
9 * Copyright 2000, 2001, 2014 Graeme W. Gill
10 * All rights reserved.
11 * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
12 * see the License.txt file for licencing details.
13 *
14 * Based on the old iccXfm class.
15 */
16
17 /*
18 * This module expands the basic icclib functionality,
19 * providing more functionality in exercising.
20 * The implementation for the three different types
21 * of profile representation, are in their own source files.
22 */
23
24 /*
25 * TTBD:
26 * Some of the error handling is crude. Shouldn't use
27 * error(), should return status.
28 *
29 */
30
31 #include <sys/types.h>
32 #include <string.h>
33 #include <ctype.h>
34 #ifdef __sun
35 #include <unistd.h>
36 #endif
37 #if defined(__IBMC__) && defined(_M_IX86)
38 #include <float.h>
39 #endif
40 #include "numlib.h"
41 #include "counters.h"
42 #include "plot.h"
43 #include "../h/sort.h"
44 #include "xicc.h" /* definitions for this library */
45
46 #define USE_CAM /* Use CIECAM02 for clipping and gamut mapping, else use Lab */
47
48 #undef DEBUG /* Plot 1d Luts */
49
50 #ifdef DEBUG
51 #include "plot.h"
52 #endif
53
54 #define MAX_INVSOLN 4
55
56 static void xicc_del(xicc *p);
57 icxLuBase * xicc_get_luobj(xicc *p, int flags, icmLookupFunc func, icRenderingIntent intent,
58 icColorSpaceSignature pcsor, icmLookupOrder order,
59 icxViewCond *vc, icxInk *ink);
60 static icxLuBase *xicc_set_luobj(xicc *p, icmLookupFunc func, icRenderingIntent intent,
61 icmLookupOrder order, int flags, int no, int nobw, cow *points,
62 icxMatrixModel *skm,
63 double dispLuminance, double wpscale,
64 // double *bpo,
65 double smooth, double avgdev,
66 double demph, icxViewCond *vc, icxInk *ink, xcal *cal, int quality);
67 static void icxLutSpaces(icxLuBase *p, icColorSpaceSignature *ins, int *inn,
68 icColorSpaceSignature *outs, int *outn,
69 icColorSpaceSignature *pcs);
70 static void icxLuSpaces(icxLuBase *p, icColorSpaceSignature *ins, int *inn,
71 icColorSpaceSignature *outs, int *outn,
72 icmLuAlgType *alg, icRenderingIntent *intt,
73 icmLookupFunc *fnc, icColorSpaceSignature *pcs);
74 static void icxLu_get_native_ranges (icxLuBase *p,
75 double *inmin, double *inmax, double *outmin, double *outmax);
76 static void icxLu_get_ranges (icxLuBase *p,
77 double *inmin, double *inmax, double *outmin, double *outmax);
78 static void icxLuEfv_wh_bk_points(icxLuBase *p, double *wht, double *blk, double *kblk);
79 int xicc_get_viewcond(xicc *p, icxViewCond *vc);
80
81 /* The different profile types are in their own source filesm */
82 /* and are included to keep their functions private. (static) */
83 #include "xmono.c"
84 #include "xmatrix.c"
85 #include "xlut.c" /* New xfit3 in & out optimising based profiles */
86 //#include "xlut1.c" /* Old xfit1 device curve based profiles */
87
88 #ifdef NT /* You'd think there might be some standards.... */
89 # ifndef __BORLANDC__
90 # define stricmp _stricmp
91 # endif
92 #else
93 # define stricmp strcasecmp
94 #endif
95 /* Utilities */
96
97 /* Return a string description of the given enumeration value */
icx2str(icmEnumType etype,int enumval)98 const char *icx2str(icmEnumType etype, int enumval) {
99
100 if (etype == icmColorSpaceSignature) {
101 if (((icColorSpaceSignature)enumval) == icxSigJabData)
102 return "Jab";
103 else if (((icColorSpaceSignature)enumval) == icxSigJChData)
104 return "JCh";
105 else if (((icColorSpaceSignature)enumval) == icxSigLChData)
106 return "LCh";
107 } else if (etype == icmRenderingIntent) {
108 if (((icRenderingIntent)enumval) == icxAppearance)
109 return "icxAppearance";
110 else if (((icRenderingIntent)enumval) == icxAbsAppearance)
111 return "icxAbsAppearance";
112 else if (((icRenderingIntent)enumval) == icxPerceptualAppearance)
113 return "icxPerceptualAppearance";
114 else if (((icRenderingIntent)enumval) == icxAbsPerceptualAppearance)
115 return "icxAbsPerceptualAppearance";
116 else if (((icRenderingIntent)enumval) == icxSaturationAppearance)
117 return "icxSaturationAppearance";
118 else if (((icRenderingIntent)enumval) == icxAbsSaturationAppearance)
119 return "icxAbsSaturationAppearance";
120 }
121 return icm2str(etype, enumval);
122 }
123
124 /* Common xicc stuff */
125
126 /* Return information about the native lut in/out colorspaces. */
127 /* Any pointer may be NULL if value is not to be returned */
128 static void
icxLutSpaces(icxLuBase * p,icColorSpaceSignature * ins,int * inn,icColorSpaceSignature * outs,int * outn,icColorSpaceSignature * pcs)129 icxLutSpaces(
130 icxLuBase *p, /* This */
131 icColorSpaceSignature *ins, /* Return input color space */
132 int *inn, /* Return number of input components */
133 icColorSpaceSignature *outs, /* Return output color space */
134 int *outn, /* Return number of output components */
135 icColorSpaceSignature *pcs /* Return PCS color space */
136 ) {
137 p->plu->lutspaces(p->plu, ins, inn, outs, outn, pcs);
138 }
139
140 /* Return information about the overall lookup in/out colorspaces, */
141 /* including allowance for any PCS override. */
142 /* Any pointer may be NULL if value is not to be returned */
143 static void
icxLuSpaces(icxLuBase * p,icColorSpaceSignature * ins,int * inn,icColorSpaceSignature * outs,int * outn,icmLuAlgType * alg,icRenderingIntent * intt,icmLookupFunc * fnc,icColorSpaceSignature * pcs)144 icxLuSpaces(
145 icxLuBase *p, /* This */
146 icColorSpaceSignature *ins, /* Return input color space */
147 int *inn, /* Return number of input components */
148 icColorSpaceSignature *outs, /* Return output color space */
149 int *outn, /* Return number of output components */
150 icmLuAlgType *alg, /* Return type of lookup algorithm used */
151 icRenderingIntent *intt, /* Return the intent implemented */
152 icmLookupFunc *fnc, /* Return the profile function being implemented */
153 icColorSpaceSignature *pcs /* Return the effective PCS */
154 ) {
155 icmLookupFunc function;
156 icColorSpaceSignature npcs; /* Native PCS */
157
158 p->plu->spaces(p->plu, NULL, inn, NULL, outn, alg, NULL, &function, &npcs, NULL);
159
160 if (intt != NULL)
161 *intt = p->intent;
162
163 if (fnc != NULL)
164 *fnc = function;
165
166 if (ins != NULL)
167 *ins = p->ins;
168
169 if (outs != NULL)
170 *outs = p->outs;
171
172 if (pcs != NULL)
173 *pcs = p->pcs;
174 }
175
176 /* Return the native (internaly visible) colorspace value ranges */
177 static void
icxLu_get_native_ranges(icxLuBase * p,double * inmin,double * inmax,double * outmin,double * outmax)178 icxLu_get_native_ranges (
179 icxLuBase *p,
180 double *inmin, double *inmax, /* Return maximum range of inspace values */
181 double *outmin, double *outmax /* Return maximum range of outspace values */
182 ) {
183 int i;
184 if (inmin != NULL) {
185 for (i = 0; i < p->inputChan; i++)
186 inmin[i] = p->ninmin[i];
187 }
188 if (inmax != NULL) {
189 for (i = 0; i < p->inputChan; i++)
190 inmax[i] = p->ninmax[i];
191 }
192 if (outmin != NULL) {
193 for (i = 0; i < p->outputChan; i++)
194 outmin[i] = p->noutmin[i];
195 }
196 if (outmax != NULL) {
197 for (i = 0; i < p->outputChan; i++)
198 outmax[i] = p->noutmax[i];
199 }
200 }
201
202 /* Return the effective (externaly visible) colorspace value ranges */
203 static void
icxLu_get_ranges(icxLuBase * p,double * inmin,double * inmax,double * outmin,double * outmax)204 icxLu_get_ranges (
205 icxLuBase *p,
206 double *inmin, double *inmax, /* Return maximum range of inspace values */
207 double *outmin, double *outmax /* Return maximum range of outspace values */
208 ) {
209 int i;
210 if (inmin != NULL) {
211 for (i = 0; i < p->inputChan; i++)
212 inmin[i] = p->inmin[i];
213 }
214 if (inmax != NULL) {
215 for (i = 0; i < p->inputChan; i++)
216 inmax[i] = p->inmax[i];
217 }
218 if (outmin != NULL) {
219 for (i = 0; i < p->outputChan; i++)
220 outmin[i] = p->outmin[i];
221 }
222 if (outmax != NULL) {
223 for (i = 0; i < p->outputChan; i++)
224 outmax[i] = p->outmax[i];
225 }
226 }
227
228 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
229 /* Routine to figure out a suitable black point for CMYK */
230
231 /* Structure to hold optimisation information */
232 typedef struct {
233 icmLuBase *p;
234 int kch; /* K channel, -1 if none */
235 double tlimit, klimit; /* Ink limit values */
236 int inn; /* Number of input channels */
237 icColorSpaceSignature outs; /* Output space */
238 double p1[3]; /* white pivot point in abs Lab */
239 double p2[3]; /* Point on vector towards black in abs Lab */
240 double toll; /* Tollerance of black direction */
241 } bpfind;
242
243 /* Optimise device values to minimise L, while remaining */
244 /* within the ink limit, and staying in line between p1 (white) and p2 (black dir) */
bpfindfunc(void * adata,double pv[])245 static double bpfindfunc(void *adata, double pv[]) {
246 bpfind *b = (bpfind *)adata;
247 double rv = 0.0;
248 double Lab[3];
249 double lr, ta, tb, terr; /* L ratio, target a, target b, target error */
250 double ovr = 0.0;
251 int e;
252
253 /* Compute amount outside total limit */
254 if (b->tlimit >= 0.0) {
255 double sum;
256 for (sum = 0.0, e = 0; e < b->inn; e++)
257 sum += pv[e];
258 if (sum > b->tlimit) {
259 ovr = sum - b->tlimit;
260 #ifdef DEBUG
261 printf("~1 total ink ovr = %f\n",ovr);
262 #endif
263 }
264 }
265
266 /* Compute amount outside black limit */
267 if (b->klimit >= 0.0 && b->kch >= 0) {
268 double kval = pv[b->kch] - b->klimit;
269 if (kval > ovr) {
270 ovr = kval;
271 #ifdef DEBUG
272 printf("~1 black ink ovr = %f\n",ovr);
273 #endif
274 }
275 }
276 /* Compute amount outside device value limits 0.0 - 1.0 */
277 {
278 double dval;
279 for (dval = -1.0, e = 0; e < b->inn; e++) {
280 if (pv[e] < 0.0) {
281 if (-pv[e] > dval)
282 dval = -pv[e];
283 } else if (pv[e] > 1.0) {
284 if ((pv[e] - 1.0) > dval)
285 dval = pv[e] - 1.0;
286 }
287 }
288 if (dval > ovr)
289 ovr = dval;
290 }
291
292 /* Compute the Lab value: */
293 b->p->lookup(b->p, Lab, pv);
294 if (b->outs == icSigXYZData)
295 icmXYZ2Lab(&icmD50, Lab, Lab);
296
297 #ifdef DEBUG
298 printf("~1 p1 = %f %f %f, p2 = %f %f %f\n",b->p1[0],b->p1[1],b->p1[2],b->p2[0],b->p2[1],b->p2[2]);
299 printf("~1 device value %f %f %f %f, Lab = %f %f %f\n",pv[0],pv[1],pv[2],pv[3],Lab[0],Lab[1],Lab[2]);
300 #endif
301
302 /* Primary aim is to minimise L value */
303 rv = Lab[0];
304
305 /* See how out of line from p1 to p2 we are */
306 lr = (Lab[0] - b->p1[0])/(b->p2[0] - b->p1[0]); /* Distance towards p2 from p1 */
307 ta = lr * (b->p2[1] - b->p1[1]) + b->p1[1]; /* Target a value */
308 tb = lr * (b->p2[2] - b->p1[2]) + b->p1[2]; /* Target b value */
309
310 terr = (ta - Lab[1]) * (ta - Lab[1])
311 + (tb - Lab[2]) * (tb - Lab[2]);
312
313 if (terr < b->toll) /* Tollerance error doesn't count until it's over tollerance */
314 terr = 0.0;
315
316 #ifdef DEBUG
317 printf("~1 target error %f\n",terr);
318 #endif
319 rv += XICC_BLACK_FIND_ABERR_WEIGHT * terr; /* Make ab match more important than min. L */
320
321 #ifdef DEBUG
322 printf("~1 out of range error %f\n",ovr);
323 #endif
324 rv += 200 * ovr;
325
326 #ifdef DEBUG
327 printf("~1 black find tc ret %f\n",rv);
328 #endif
329 return rv;
330 }
331
332 /* Try and compute a real black point in XYZ given an iccLu, */
333 /* and also return the K only black or the normal black if the device doesn't have K */
334 /* black[] will be unchanged if black cannot be computed. */
335 /* Note that the black point will be in the space of the Lu */
336 /* converted to XYZ, so will have the Lu's intent etc. */
337 /* (Note that this is duplicated in xlut.c set_icxLuLut() !!!) */
icxLu_comp_bk_point(icxLuBase * x,int gblk,double * white,double * black,double * kblack)338 static void icxLu_comp_bk_point(
339 icxLuBase *x,
340 int gblk, /* If nz, compute black if possible. */
341 double *white, /* XYZ Input, used for computing black */
342 double *black, /* XYZ Input & Output. Set if gblk NZ and can be computed */
343 double *kblack /* XYZ Output. Looked up if possible or set to black[] otherwise */
344 ) {
345 icmLuBase *p = x->plu;
346 icmLuBase *op = p; /* Original icmLu, in case we replace p */
347 icc *icco = p->icp;
348 icmHeader *h = icco->header;
349 icColorSpaceSignature ins, outs;
350 int inn, outn;
351 icmLuAlgType alg;
352 icRenderingIntent intt;
353 icmLookupFunc fnc;
354 icmLookupOrder ord;
355 int kch = -1;
356 double dblack[MAX_CHAN]; /* device black value */
357 int e;
358
359 #ifdef DEBUG
360 printf("~1 icxLu_comp_bk_point() called, gblk %d, white = %s, black = %s\n",gblk,icmPdv(3, white),icmPdv(3,black));
361 #endif
362 /* Default return incoming black as K only black */
363 kblack[0] = black[0];
364 kblack[1] = black[1];
365 kblack[2] = black[2];
366
367 /* Get the effective characteristics of the Lu */
368 p->spaces(p, &ins, &inn, &outs, &outn, &alg, &intt, &fnc, NULL, &ord);
369
370 if (fnc == icmBwd) { /* Hmm. We've got PCS to device, and we want device to PCS. */
371
372 /* Strictly speaking this is a dubious approach, since for a cLut profile */
373 /* the B2A table could make the effective white and black points */
374 /* anything it likes, and they don't have to match what the corresponding */
375 /* A2B table does. In our usage it's probably OK, since we tend */
376 /* to use colorimetric B2A */
377 #ifdef DEBUG
378 printf("~1 getting icmFwd\n");
379 #endif
380 if ((p = icco->get_luobj(icco, icmFwd, intt, ins, ord)) == NULL)
381 error("icxLu_comp_bk_point: assert: getting Fwd Lookup failed!");
382
383 p->spaces(p, &ins, &inn, &outs, &outn, &alg, &intt, &fnc, NULL, &ord);
384 }
385
386 if (outs != icSigXYZData && outs != icSigLabData) {
387 error("icxLu_comp_bk_point: assert: icc Lu output is not XYZ or Lab!, outs = 0x%x, ");
388 }
389
390 #ifdef DEBUG
391 printf("~1 icxLu_comp_bk_point called for inn = %d, ins = %s\n", inn, icx2str(icmColorSpaceSignature,ins));
392 #endif
393
394 switch (ins) {
395
396 case icSigXYZData:
397 case icSigLabData:
398 case icSigLuvData:
399 case icSigYxyData:
400 #ifdef DEBUG
401 printf("~1 Assuming CIE colorspace black is 0.0\n");
402 #endif
403 if (gblk) {
404 for (e = 0; e < inn; e++)
405 black[0] = 0.0;
406 }
407 kblack[0] = black[0];
408 kblack[1] = black[1];
409 kblack[2] = black[2];
410 return;
411
412 case icSigRgbData:
413 #ifdef DEBUG
414 printf("~1 RGB:\n");
415 #endif
416 for (e = 0; e < inn; e++)
417 dblack[e] = 0.0;
418 break;
419
420 case icSigGrayData: { /* Could be additive or subtractive */
421 double dval[1];
422 double minv[3], maxv[3];
423 #ifdef DEBUG
424 printf("~1 Gray:\n");
425 #endif
426 /* Check out 0 and 100% colorant */
427 dval[0] = 0.0;
428 p->lookup(p, minv, dval);
429 if (outs == icSigXYZData)
430 icmXYZ2Lab(&icmD50, minv, minv);
431 dval[0] = 1.0;
432 p->lookup(p, maxv, dval);
433 if (outs == icSigXYZData)
434 icmXYZ2Lab(&icmD50, maxv, maxv);
435
436 if (minv[0] < maxv[0])
437 dblack[0] = 0.0;
438 else
439 dblack[0] = 1.0;
440 }
441 break;
442
443 case icSigCmyData:
444 for (e = 0; e < inn; e++)
445 dblack[e] = 1.0;
446 break;
447
448 case icSigCmykData:
449 #ifdef DEBUG
450 printf("~1 CMYK:\n");
451 #endif
452 kch = 3;
453 dblack[0] = 0.0;
454 dblack[1] = 0.0;
455 dblack[2] = 0.0;
456 dblack[3] = 1.0;
457 if (alg == icmLutType) {
458 icxLuLut *pp = (icxLuLut *)x;
459
460 if (pp->ink.tlimit >= 0.0)
461 dblack[kch] = pp->ink.tlimit;
462 };
463 break;
464
465 /* Use a heursistic. */
466 /* This duplicates code in icxGetLimits() :-( */
467 /* Colorant guessing should go in icclib ? */
468 case icSig2colorData:
469 case icSig3colorData:
470 case icSig4colorData:
471 case icSig5colorData:
472 case icSig6colorData:
473 case icSig7colorData:
474 case icSig8colorData:
475 case icSig9colorData:
476 case icSig10colorData:
477 case icSig11colorData:
478 case icSig12colorData:
479 case icSig13colorData:
480 case icSig14colorData:
481 case icSig15colorData:
482 case icSigMch5Data:
483 case icSigMch6Data:
484 case icSigMch7Data:
485 case icSigMch8Data: {
486 double dval[MAX_CHAN];
487 double ncval[3];
488 double cvals[MAX_CHAN][3];
489 int nlighter, ndarker;
490
491 /* Decide if the colorspace is additive or subtractive */
492 #ifdef DEBUG
493 printf("~1 N channel:\n");
494 #endif
495
496 /* First the no colorant value */
497 for (e = 0; e < inn; e++)
498 dval[e] = 0.0;
499 p->lookup(p, ncval, dval);
500 if (outs == icSigXYZData)
501 icmXYZ2Lab(&icmD50, ncval, ncval);
502
503 /* Then all the colorants */
504 nlighter = ndarker = 0;
505 for (e = 0; e < inn; e++) {
506 dval[e] = 1.0;
507 p->lookup(p, cvals[e], dval);
508 if (outs == icSigXYZData)
509 icmXYZ2Lab(&icmD50, cvals[e], cvals[e]);
510 dval[e] = 0.0;
511 if (fabs(cvals[e][0] - ncval[0]) > 5.0) {
512 if (cvals[e][0] > ncval[0])
513 nlighter++;
514 else
515 ndarker++;
516 }
517 }
518 if (ndarker == 0 && nlighter > 0) { /* Assume additive */
519 for (e = 0; e < inn; e++)
520 dblack[e] = 0.0;
521 #ifdef DEBUG
522 printf("~1 N channel is additive:\n");
523 #endif
524
525 } else if (ndarker > 0 && nlighter == 0) { /* Assume subtractive. */
526 double pbk[3] = { 0.0,0.0,0.0 }; /* Perfect black */
527 double smd = 1e10; /* Smallest distance */
528
529 #ifdef DEBUG
530 printf("~1 N channel is subtractive:\n");
531 #endif
532 /* See if we can guess the black channel */
533 for (e = 0; e < inn; e++) {
534 double tt;
535 tt = icmNorm33sq(pbk, cvals[e]);
536 if (tt < smd) {
537 smd = tt;
538 kch = e;
539 }
540 }
541 /* See if the black seems sane */
542 if (cvals[kch][0] > 40.0
543 || fabs(cvals[kch][1]) > 10.0
544 || fabs(cvals[kch][2]) > 10.0) {
545 if (p != op)
546 p->del(p);
547 #ifdef DEBUG
548 printf("~1 black doesn't look sanem so assume nothing\n");
549 #endif
550 return; /* Assume nothing */
551 }
552
553 /* Chosen kch as black */
554 for (e = 0; e < inn; e++)
555 dblack[e] = 0.0;
556 dblack[kch] = 1.0;
557 if (alg == icmLutType) {
558 icxLuLut *pp = (icxLuLut *)x;
559
560 if (pp->ink.tlimit >= 0.0)
561 dblack[kch] = pp->ink.tlimit;
562 };
563 #ifdef DEBUG
564 printf("~1 N channel K = chan %d\n",kch);
565 #endif
566 } else {
567 if (p != op)
568 p->del(p);
569 #ifdef DEBUG
570 printf("~1 can't figure if additive or subtractive, so assume nothing\n");
571 #endif
572 return; /* Assume nothing */
573 }
574 }
575 break;
576
577 default:
578 #ifdef DEBUG
579 printf("~1 unhandled colorspace, so assume nothing\n");
580 #endif
581 if (p != op)
582 p->del(p);
583 return; /* Don't do anything */
584 }
585
586 /* Lookup the K only value */
587 if (kch >= 0) {
588 p->lookup(p, kblack, dblack);
589
590 /* We always return XYZ */
591 if (outs == icSigLabData)
592 icmLab2XYZ(&icmD50, kblack, kblack);
593 }
594
595 if (gblk == 0) { /* That's all we have to do */
596 #ifdef DEBUG
597 printf("~1 gblk == 0, so only return kblack\n");
598 #endif
599 if (p != op)
600 p->del(p);
601 return;
602 }
603
604 /* Lookup the device black or K only value as a default */
605 p->lookup(p, black, dblack); /* May be XYZ or Lab */
606 #ifdef DEBUG
607 printf("~1 Got default lu black %f %f %f, kch = %d\n", black[0],black[1],black[2],kch);
608 #endif
609
610 /* !!! Hmm. For CMY and RGB we are simply using the device */
611 /* combination values as the black point. In reality we might */
612 /* want to have the option of using a neutral black point, */
613 /* just like CMYK ?? */
614
615 if (kch >= 0) { /* The space is subtractive with a K channel. */
616 /* If XICC_NEUTRAL_CMYK_BLACK then locate the darkest */
617 /* CMYK within limits with the same chromaticity as the white point, */
618 /* otherwise locate the device value within the ink limits that is */
619 /* in the direction of the K channel */
620 bpfind bfs; /* Callback context */
621 double sr[MXDO]; /* search radius */
622 double tt[MXDO]; /* Temporary */
623 double rs0[MXDO], rs1[MXDO]; /* Random start candidates */
624 int trial;
625 double brv;
626
627 /* Setup callback function context */
628 bfs.p = p;
629 bfs.inn = inn;
630 bfs.outs = outs;
631
632 bfs.kch = kch;
633 bfs.tlimit = -1.0;
634 bfs.klimit = -1.0;
635 bfs.toll = XICC_BLACK_POINT_TOLL;
636
637 if (alg == icmLutType) {
638 icxLuLut *pp = (icxLuLut *)x;
639
640 pp->kch = kch;
641 bfs.tlimit = pp->ink.tlimit;
642 bfs.klimit = pp->ink.klimit;
643 #ifdef DEBUG
644 printf("~1 tlimit = %f, klimit = %f\n",bfs.tlimit,bfs.klimit);
645 #endif
646 };
647
648 #ifdef XICC_NEUTRAL_CMYK_BLACK
649 #ifdef DEBUG
650 printf("~1 Searching for neutral black\n");
651 #endif
652 /* white has been given to us in XYZ */
653 icmXYZ2Lab(&icmD50, bfs.p1, white); /* pivot Lab */
654 icmCpy3(bfs.p2, white); /* temp white XYZ */
655 icmScale3(bfs.p2, bfs.p2, 0.02); /* Scale white XYZ towards 0,0,0 */
656 icmXYZ2Lab(&icmD50, bfs.p2, bfs.p2); /* Convert black direction to Lab */
657
658 #else /* Use K directin black */
659 #ifdef DEBUG
660 printf("~1 Searching for K direction black\n");
661 #endif
662 icmXYZ2Lab(&icmD50, bfs.p1, white); /* Pivot */
663
664 /* Now figure abs Lab value of K only, as the direction */
665 /* to use for the rich black. */
666 for (e = 0; e < inn; e++)
667 dblack[e] = 0.0;
668 if (bfs.klimit < 0.0)
669 dblack[kch] = 1.0;
670 else
671 dblack[kch] = bfs.klimit; /* K value */
672
673 p->lookup(p, black, dblack);
674
675 if (outs == icSigXYZData) {
676 icmXYZ2Lab(&icmD50, bfs.p2, black); /* K direction */
677 } else {
678 icmAry2Ary(bfs.p2, black);
679 }
680 #endif
681
682 #ifdef DEBUG
683 printf("~1 Lab pivot %f %f %f, Lab K direction %f %f %f\n",bfs.p1[0],bfs.p1[1],bfs.p1[2],bfs.p2[0],bfs.p2[1],bfs.p2[2]);
684 #endif
685 /* Set the random start 0 location as 000K */
686 /* and the random start 1 location as CMY0 */
687 {
688 double tt;
689
690 for (e = 0; e < inn; e++)
691 dblack[e] = rs0[e] = 0.0;
692 if (bfs.klimit < 0.0)
693 dblack[kch] = rs0[kch] = 1.0;
694 else
695 dblack[kch] = rs0[kch] = bfs.klimit; /* K value */
696
697 if (bfs.tlimit < 0.0)
698 tt = 1.0;
699 else
700 tt = bfs.tlimit/(inn - 1.0);
701 for (e = 0; e < inn; e++)
702 rs1[e] = tt;
703 rs1[kch] = 0.0; /* K value */
704 }
705
706 /* Start with the K only as the current best value */
707 brv = bpfindfunc((void *)&bfs, dblack);
708 #ifdef DEBUG
709 printf("~1 initial brv for K only = %f\n",brv);
710 #endif
711
712 /* Find the device black point using optimization */
713 /* Do several trials to avoid local minima. */
714 rand32(0x12345678); /* Make trial values deterministic */
715 for (trial = 0; trial < 200; trial++) {
716 double rv; /* Temporary */
717
718 /* Start first trial at 000K */
719 if (trial == 0) {
720 for (e = 0; e < inn; e++) {
721 tt[e] = rs0[e];
722 sr[e] = 0.1;
723 }
724
725 } else {
726 /* Base is random between 000K and CMY0: */
727 if (trial < 100) {
728 rv = d_rand(0.0, 1.0);
729 for (e = 0; e < inn; e++) {
730 tt[e] = rv * rs0[e] + (1.0 - rv) * rs1[e];
731 sr[e] = 0.1;
732 }
733 /* Base on current best */
734 } else {
735 for (e = 0; e < inn; e++) {
736 tt[e] = dblack[e];
737 sr[e] = 0.1;
738 }
739 }
740
741 /* Then add random start offset */
742 for (rv = 0.0, e = 0; e < inn; e++) {
743 tt[e] += d_rand(-0.5, 0.5);
744 if (tt[e] < 0.0)
745 tt[e] = 0.0;
746 else if (tt[e] > 1.0)
747 tt[e] = 1.0;
748 }
749 }
750
751 /* Clip black */
752 if (bfs.klimit >= 0.0 && tt[kch] > bfs.klimit)
753 tt[kch] = bfs.klimit;
754
755 /* Compute amount outside total limit */
756 if (bfs.tlimit >= 0.0) {
757 for (rv = 0.0, e = 0; e < inn; e++)
758 rv += tt[e];
759
760 if (rv > bfs.tlimit) {
761 rv /= (double)inn;
762 for (e = 0; e < inn; e++)
763 tt[e] -= rv;
764 }
765 }
766
767 if (powell(&rv, inn, tt, sr, 0.000001, 1000, bpfindfunc,
768 (void *)&bfs, NULL, NULL) == 0) {
769 #ifdef DEBUG
770 printf("~1 trial %d, rv %f bp %f %f %f %f\n",trial,rv,tt[0],tt[1],tt[2],tt[3]);
771 #endif
772 if (rv < brv) {
773 #ifdef DEBUG
774 printf("~1 new best\n");
775 #endif
776 brv = rv;
777 for (e = 0; e < inn; e++)
778 dblack[e] = tt[e];
779 }
780 }
781 }
782 if (brv > 1000.0)
783 error("icxLu_comp_bk_point: Black point powell failed");
784
785 for (e = 0; e < inn; e++) { /* Make sure device values are in range */
786 if (dblack[e] < 0.0)
787 dblack[e] = 0.0;
788 else if (dblack[e] > 1.0)
789 dblack[e] = 1.0;
790 }
791 /* Now have device black in dblack[] */
792 #ifdef DEBUG
793 printf("~1 got device black %f %f %f %f\n",dblack[0], dblack[1], dblack[2], dblack[3]);
794 #endif
795
796 p->lookup(p, black, dblack); /* Convert to PCS */
797 }
798
799 if (p != op)
800 p->del(p);
801
802 /* We always return XYZ */
803 if (outs == icSigLabData)
804 icmLab2XYZ(&icmD50, black, black);
805
806 #ifdef DEBUG
807 printf("~1 returning %f %f %f\n", black[0], black[1], black[2]);
808 #endif
809
810 return;
811 }
812
813 /* - - - - - - - - - - - - - - - - - - - - - - - */
814
815 /* Return the media white and black points */
816 /* in the xlu effective PCS colorspace. Pointers may be NULL. */
817 /* (ie. these will be relative values for relative intent etc.) */
icxLuEfv_wh_bk_points(icxLuBase * p,double * wht,double * blk,double * kblk)818 static void icxLuEfv_wh_bk_points(
819 icxLuBase *p,
820 double *wht,
821 double *blk,
822 double *kblk /* K only black */
823 ) {
824 double white[3], black[3], kblack[3];
825
826 /* Get the Lu PCS converted to XYZ icc black and white points in XYZ */
827 if (p->plu->lu_wh_bk_points(p->plu, white, black)) {
828 /* Black point is assumed. We should determine one instead. */
829 /* Lookup K only black too */
830 icxLu_comp_bk_point(p, 1, white, black, kblack);
831
832 } else {
833 /* Lookup a possible K only black */
834 icxLu_comp_bk_point(p, 0, white, black, kblack);
835 }
836
837 //printf("~1 white %f %f %f, black %f %f %f, kblack %f %f %f\n",white[0],white[1],white[2],black[0],black[1],black[2],kblack[0],kblack[1],kblack[2]);
838 /* Convert to possibl xicc override PCS */
839 switch (p->pcs) {
840 case icSigXYZData:
841 break; /* Don't have to do anyting */
842 case icSigLabData:
843 icmXYZ2Lab(&icmD50, white, white); /* Convert from XYZ to Lab */
844 icmXYZ2Lab(&icmD50, black, black);
845 icmXYZ2Lab(&icmD50, kblack, kblack);
846 break;
847 case icxSigJabData:
848 p->cam->XYZ_to_cam(p->cam, white, white); /* Convert from XYZ to Jab */
849 p->cam->XYZ_to_cam(p->cam, black, black);
850 p->cam->XYZ_to_cam(p->cam, kblack, kblack);
851 break;
852 default:
853 break;
854 }
855
856 //printf("~1 icxLuEfv_wh_bk_points: pcsor %s White %f %f %f, Black %f %f %f\n", icx2str(icmColorSpaceSignature,p->pcs), white[0], white[1], white[2], black[0], black[1], black[2]);
857 if (wht != NULL) {
858 wht[0] = white[0];
859 wht[1] = white[1];
860 wht[2] = white[2];
861 }
862
863 if (blk != NULL) {
864 blk[0] = black[0];
865 blk[1] = black[1];
866 blk[2] = black[2];
867 }
868
869 if (kblk != NULL) {
870 kblk[0] = kblack[0];
871 kblk[1] = kblack[1];
872 kblk[2] = kblack[2];
873 }
874 }
875
876 /* Create an instance of an xicc object */
new_xicc(icc * picc)877 xicc *new_xicc(
878 icc *picc /* icc we are expanding */
879 ) {
880 xicc *p;
881 if ((p = (xicc *) calloc(1,sizeof(xicc))) == NULL)
882 return NULL;
883 p->pp = picc;
884 p->del = xicc_del;
885 p->get_luobj = xicc_get_luobj;
886 p->set_luobj = xicc_set_luobj;
887 p->get_viewcond = xicc_get_viewcond;
888
889 /* Create an xcal if there is the right tag in the profile */
890 p->cal = xiccReadCalTag(p->pp);
891 p->nodel_cal = 0; /* We created it, we will delete it */
892
893 return p;
894 }
895
896 /* Do away with the xicc (but not the icc!) */
xicc_del(xicc * p)897 static void xicc_del(
898 xicc *p
899 ) {
900 if (p->cal != NULL && p->nodel_cal == 0)
901 p->cal->del(p->cal);
902 free (p);
903 }
904
905 /* return nz if the intent implies Jab space */
xiccIsIntentJab(icRenderingIntent intent)906 int xiccIsIntentJab(icRenderingIntent intent) {
907
908 if (intent == icxAppearance
909 || intent == icxAbsAppearance
910 || intent == icxPerceptualAppearance
911 || intent == icxAbsPerceptualAppearance
912 || intent == icxSaturationAppearance
913 || intent == icxAbsSaturationAppearance)
914 return 1;
915 return 0;
916 }
917
918 /* Return an expanded lookup object, initialised */
919 /* from the icc. */
920 /* Return NULL on error, check errc+err for reason. */
921 /* Set the pcsor & intent to consistent and values if */
922 /* Jab and/or icxAppearance has been requested. */
923 /* Create the underlying icm lookup object that is used */
924 /* to create and implement the icx one. The icm will be used */
925 /* to translate from native to effective PCS, unless the */
926 /* effective PCS is Jab, in which case the icm will be set to */
927 /* have an effective PCS of XYZ. Since native<->effecive PCS conversion */
928 /* is done at the to/from_abs() stage, none of this affects the individual */
929 /* conversion steps, which will all talk the native PCS (unless merged). */
xicc_get_luobj(xicc * p,int flags,icmLookupFunc func,icRenderingIntent intent,icColorSpaceSignature pcsor,icmLookupOrder order,icxViewCond * vc,icxInk * ink)930 icxLuBase *xicc_get_luobj(
931 xicc *p, /* this */
932 int flags, /* clip, merge flags */
933 icmLookupFunc func, /* Functionality */
934 icRenderingIntent intent, /* Intent */
935 icColorSpaceSignature pcsor,/* PCS override (0 = def) */
936 icmLookupOrder order, /* Search Order */
937 icxViewCond *vc, /* Viewing Condition (may be NULL if pcsor is not CIECAM) */
938 icxInk *ink /* inking details (NULL for default) */
939 ) {
940 icmLuBase *plu;
941 icxLuBase *xplu;
942 icmLuAlgType alg;
943 icRenderingIntent n_intent = intent; /* Native Intent to request */
944 icColorSpaceSignature n_pcs = icmSigDefaultData; /* Native PCS to request */
945
946 //printf("~1 xicc_get_luobj got intent '%s' and pcsor '%s'\n",icx2str(icmRenderingIntent,intent),icx2str(icmColorSpaceSignature,pcsor));
947
948 /* Ensure that appropriate PCS is slected for an appearance intent */
949 if (xiccIsIntentJab(intent)) {
950 pcsor = icxSigJabData;
951 //printf("~1 pcsor = %s\n",tag2str(pcsor));
952
953 /* Translate non-Jab intents to the equivalent appearance "intent" if pcsor == Jab. */
954 /* This is how we get these when the UI's don't list all the apperances intents, */
955 /* we select the analogous non-apperance intent with pcsor = Jab. */
956 /* Note that Abs/non-abs selects between Apperance and AbsAppearance. */
957 } else if (pcsor == icxSigJabData) {
958 if (intent == icRelativeColorimetric)
959 intent = icxAppearance;
960 else if (intent == icAbsoluteColorimetric)
961 intent = icxAbsAppearance;
962 else if (intent == icPerceptual)
963 intent = icxPerceptualAppearance;
964 else if (intent == icmAbsolutePerceptual)
965 intent = icxAbsPerceptualAppearance;
966 else if (intent == icSaturation)
967 intent = icxSaturationAppearance;
968 else if (intent == icmAbsoluteSaturation)
969 intent = icxAbsSaturationAppearance;
970 else
971 intent = icxAppearance;
972 }
973 //printf("~1 intent = %s\n",tag2str(intent));
974
975 /* Translate intent asked for into intent needed in icclib */
976 if (intent == icxAppearance
977 || intent == icxAbsAppearance)
978 n_intent = icAbsoluteColorimetric;
979 else if (intent == icxPerceptualAppearance
980 || intent == icxAbsPerceptualAppearance)
981 n_intent = icmAbsolutePerceptual;
982 else if (intent == icxSaturationAppearance
983 || intent == icxAbsSaturationAppearance)
984 n_intent = icmAbsoluteSaturation;
985 //printf("~1 n_intent = %s\n",tag2str(n_intent));
986
987 if (pcsor != icmSigDefaultData)
988 n_pcs = pcsor; /* There is an icclib override */
989
990 if (pcsor == icxSigJabData) /* xicc override */
991 n_pcs = icSigXYZData; /* Translate to XYZ */
992
993 //printf("~1 xicc_get_luobj processed intent %s and pcsor %s\n",icx2str(icmRenderingIntent,intent),icx2str(icmColorSpaceSignature,pcsor));
994 //printf("~1 xicc_get_luobj icclib intent %s and pcsor %s\n",icx2str(icmRenderingIntent,n_intent),icx2str(icmColorSpaceSignature,n_pcs));
995 /* Get icclib lookup object */
996 if ((plu = p->pp->get_luobj(p->pp, func, n_intent, n_pcs, order)) == NULL) {
997 p->errc = p->pp->errc; /* Copy error */
998 strcpy(p->err, p->pp->err);
999 return NULL;
1000 }
1001
1002 /* Figure out what the algorithm is */
1003 plu->spaces(plu, NULL, NULL, NULL, NULL, &alg, NULL, NULL, &n_pcs, NULL);
1004
1005 /* make sure its "Abs CAM" */
1006 if (vc!= NULL
1007 && (intent == icxAbsAppearance
1008 || intent == icxAbsPerceptualAppearance
1009 || intent == icxAbsSaturationAppearance)) { /* make sure its "Abs CAM" */
1010 //printf("~1 xicc_get_luobj using absolute apperance space with white = D50\n");
1011 /* Set white point and flare color to D50 */
1012 /* (Hmm. This doesn't match what happens within collink with absolute intent!!) */
1013 vc->Wxyz[0] = icmD50.X/icmD50.Y;
1014 vc->Wxyz[1] = icmD50.Y/icmD50.Y; // Normalise white reference to Y = 1 ?
1015 vc->Wxyz[2] = icmD50.Z/icmD50.Y;
1016
1017 vc->Gxyz[0] = icmD50.X;
1018 vc->Gxyz[1] = icmD50.Y;
1019 vc->Gxyz[2] = icmD50.Z;
1020 }
1021
1022 /* Call xiccLu wrapper creation */
1023 switch (alg) {
1024 case icmMonoFwdType:
1025 xplu = new_icxLuMono(p, flags, plu, func, intent, pcsor, vc, 0);
1026 break;
1027 case icmMonoBwdType:
1028 xplu = new_icxLuMono(p, flags, plu, func, intent, pcsor, vc, 1);
1029 break;
1030 case icmMatrixFwdType:
1031 xplu = new_icxLuMatrix(p, flags, plu, func, intent, pcsor, vc, 0);
1032 break;
1033 case icmMatrixBwdType:
1034 xplu = new_icxLuMatrix(p, flags, plu, func, intent, pcsor, vc, 1);
1035 break;
1036 case icmLutType:
1037 xplu = new_icxLuLut(p, flags, plu, func, intent, pcsor, vc, ink);
1038 break;
1039 default:
1040 xplu = NULL;
1041 break;
1042 }
1043
1044 return xplu;
1045 }
1046
1047
1048 /* Return an expanded lookup object, initialised */
1049 /* from the icc, and then overwritten by a conversion */
1050 /* created from the supplied scattered data points. */
1051 /* The Lut is assumed to be a device -> native PCS profile. */
1052 /* If the SET_WHITE and/or SET_BLACK flags are set, */
1053 /* discover the white/black point, set it in the icc, */
1054 /* and make the Lut relative to them. */
1055 /* Return NULL on error, check errc+err for reason */
xicc_set_luobj(xicc * p,icmLookupFunc func,icRenderingIntent intent,icmLookupOrder order,int flags,int no,int nobw,cow * points,icxMatrixModel * skm,double dispLuminance,double wpscale,double smooth,double avgdev,double demph,icxViewCond * vc,icxInk * ink,xcal * cal,int quality)1056 static icxLuBase *xicc_set_luobj(
1057 xicc *p, /* this */
1058 icmLookupFunc func, /* Functionality */
1059 icRenderingIntent intent, /* Intent */
1060 icmLookupOrder order, /* Search Order */
1061 int flags, /* white/black point, verbose flags etc. */
1062 int no, /* Number of points */
1063 int nobw, /* Number of points to look for white & black patches in */
1064 cow *points, /* Array of input points in target PCS space */
1065 icxMatrixModel *skm, /* Optional skeleton model (used for input profiles) */
1066 double dispLuminance, /* > 0.0 if display luminance value and is known */
1067 double wpscale, /* > 0.0 if input white point is to be scaled */
1068 //double *bpo, /* != NULL for black point override XYZ */
1069 double smooth, /* RSPL smoothing factor, -ve if raw */
1070 double avgdev, /* reading Average Deviation as a proportion of the input range */
1071 double demph, /* dark emphasis factor for cLUT grid res. */
1072 icxViewCond *vc, /* Viewing Condition (NULL if not using CAM) */
1073 icxInk *ink, /* inking details (NULL for default) */
1074 xcal *cal, /* Optional cal, will override any existing (not deleted with xicc)*/
1075 int quality /* Quality metric, 0..3 */
1076 ) {
1077 icmLuBase *plu;
1078 icxLuBase *xplu = NULL;
1079 icmLuAlgType alg;
1080
1081 if (cal != NULL) {
1082 if (p->cal != NULL && p->nodel_cal == 0)
1083 p->cal->del(p->cal);
1084 p->cal = cal;
1085 p->nodel_cal = 1; /* We were given it, so don't delete it */
1086 }
1087
1088 if (func != icmFwd) {
1089 p->errc = 1;
1090 sprintf(p->err,"Can only create Device->PCS profiles from scattered data.");
1091 xplu = NULL;
1092 return xplu;
1093 }
1094
1095 /* Get icclib lookup object */
1096 if ((plu = p->pp->get_luobj(p->pp, func, intent, 0, order)) == NULL) {
1097 p->errc = p->pp->errc; /* Copy error */
1098 strcpy(p->err, p->pp->err);
1099 return NULL;
1100 }
1101
1102 /* Figure out what the algorithm is */
1103 plu->spaces(plu, NULL, NULL, NULL, NULL, &alg, NULL, NULL, NULL, NULL);
1104
1105 /* Call xiccLu wrapper creation */
1106 switch (alg) {
1107 case icmMonoFwdType:
1108 p->errc = 1;
1109 sprintf(p->err,"Setting Monochrome Fwd profile from scattered data not supported.");
1110 plu->del(plu);
1111 xplu = NULL; /* Not supported yet */
1112 break;
1113
1114 case icmMatrixFwdType:
1115 if (smooth < 0.0)
1116 smooth = -smooth;
1117 xplu = set_icxLuMatrix(p, plu, flags, no, nobw, points, skm, dispLuminance, wpscale,
1118 // bpo,
1119 quality, smooth);
1120 break;
1121
1122 case icmLutType:
1123 /* ~~~ Should add check that it is a fwd profile ~~~ */
1124 xplu = set_icxLuLut(p, plu, func, intent, flags, no, nobw, points, skm, dispLuminance,
1125 wpscale,
1126 // bpo,
1127 smooth, avgdev, demph, vc, ink, quality);
1128 break;
1129
1130 default:
1131 break;
1132 }
1133
1134 return xplu;
1135 }
1136
1137 /* ------------------------------------------------------ */
1138 /* Viewing Condition Parameter stuff */
1139
1140 #ifdef NEVER /* Not currently used */
1141
1142 /* Guess viewing parameters from the technology signature */
guess_from_techsig(icTechnologySignature tsig,double * Ybp)1143 static void guess_from_techsig(
1144 icTechnologySignature tsig,
1145 double *Ybp
1146 ) {
1147 double Yb = -1.0;
1148
1149 switch (tsig) {
1150 /* These are all inputing either a representation of */
1151 /* a natural scene captured on another medium, or are assuming */
1152 /* that the medium is the original. A _good_ system would */
1153 /* let the user indicate which is the case. */
1154 case icSigReflectiveScanner:
1155 case icSigFilmScanner:
1156 Yb = 0.2;
1157 break;
1158
1159 /* Direct scene to value devices. */
1160 case icSigDigitalCamera:
1161 case icSigVideoCamera:
1162 Yb = 0.2;
1163 break;
1164
1165 /* Emmisive displays. */
1166 /* We could try tweaking the white point on the assumption */
1167 /* that the viewer will be adapted to a combination of both */
1168 /* the CRT white point, and the ambient light. */
1169 case icSigVideoMonitor:
1170 case icSigCRTDisplay:
1171 case icSigPMDisplay:
1172 case icSigAMDisplay:
1173 Yb = 0.2;
1174 break;
1175
1176 /* Photo CD has its own viewing definitions */
1177 /* (It represents original scene colors) */
1178 case icSigPhotoCD:
1179 Yb = 0.2;
1180 break;
1181
1182 /* Projection devices, either direct, or */
1183 /* via another intermediate medium. */
1184 case icSigProjectionTelevision:
1185 Yb = 0.1; /* Assume darkened room, little background */
1186 break;
1187 case icSigFilmWriter:
1188 Yb = 0.0; /* Assume a dark room - no background */
1189 break;
1190
1191 /* Printed media devices. */
1192 case icSigInkJetPrinter:
1193 case icSigThermalWaxPrinter:
1194 case icSigElectrophotographicPrinter:
1195 case icSigElectrostaticPrinter:
1196 case icSigDyeSublimationPrinter:
1197 case icSigPhotographicPaperPrinter:
1198 case icSigPhotoImageSetter:
1199 case icSigGravure:
1200 case icSigOffsetLithography:
1201 case icSigSilkscreen:
1202 case icSigFlexography:
1203 Yb = 0.2;
1204 break;
1205
1206 default:
1207 Yb = 0.2;
1208 }
1209
1210 if (Ybp != NULL)
1211 *Ybp = Yb;
1212 }
1213
1214 #endif /* NEVER */
1215
1216
1217 /* See if we can read or guess the viewing conditions */
1218 /* for an ICC profile. */
1219 /* Return value 0 if it is well defined */
1220 /* Return value 1 if it is a guess */
1221 /* Return value 2 if it is not possible/appropriate */
xicc_get_viewcond(xicc * p,icxViewCond * vc)1222 int xicc_get_viewcond(
1223 xicc *p, /* Expanded profile we're working with */
1224 icxViewCond *vc /* Viewing parameters to return */
1225 ) {
1226 icc *pp = p->pp; /* Base ICC */
1227
1228 /* Numbers we're trying to find */
1229 ViewingCondition Ev = vc_none;
1230 double Wxyz[3] = {-1.0, -1.0, -1.0}; /* Adapting white color */
1231 double La = -1.0; /* Adapting/Surround luminance */
1232 double Ixyz[3] = {-1.0, -1.0, -1.0}; /* Illuminant color */
1233 double Li = -1.0; /* Illuminant luminance */
1234 double Lb = -1.0; /* Backgrount luminance */
1235 double Yb = -1.0; /* Background relative luminance to Lv */
1236 double Lve = -1.0; /* Emissive device image luminance */
1237 double Lvr = -1.0; /* Reflective device image luminance */
1238 double Lv = -1.0; /* Device image luminance */
1239 double Yf = -1.0; /* Flare relative luminance to Lv */
1240 double Yg = -1.0; /* Glare relative luminance to La */
1241 double Gxyz[3] = {-1.0, -1.0, -1.0}; /* Glare color */
1242 icTechnologySignature tsig = icMaxEnumTechnology; /* Technology Signature */
1243 icProfileClassSignature devc = icMaxEnumClass;
1244 int trans = -1; /* Set to 0 if not transparency, 1 if it is */
1245
1246 /* Collect all the information we can find */
1247
1248 /* Emmisive devices image white luminance */
1249 {
1250 icmXYZArray *luminanceTag;
1251
1252 if ((luminanceTag = (icmXYZArray *)pp->read_tag(pp, icSigLuminanceTag)) != NULL
1253 && luminanceTag->ttype == icSigXYZType && luminanceTag->size >= 1) {
1254 Lve = luminanceTag->data[0].Y; /* Copy structure */
1255 }
1256 }
1257
1258 /* Flare: */
1259 {
1260 icmMeasurement *ro;
1261
1262 if ((ro = (icmMeasurement *)pp->read_tag(pp, icSigMeasurementTag)) != NULL
1263 && ro->ttype == icSigMeasurementType) {
1264
1265 Yf = 0.0 * ro->flare; // ?????
1266 Yg = 1.0 * ro->flare; // ?????
1267 /* ro->illuminant ie D50, D65, D93, A etc. */
1268 }
1269 }
1270
1271 /* Media White Point */
1272 {
1273 icmXYZArray *whitePointTag;
1274
1275 if ((whitePointTag = (icmXYZArray *)pp->read_tag(pp, icSigMediaWhitePointTag)) != NULL
1276 && whitePointTag->ttype == icSigXYZType && whitePointTag->size >= 1) {
1277 Wxyz[0] = whitePointTag->data[0].X;
1278 Wxyz[1] = whitePointTag->data[0].Y;
1279 Wxyz[2] = whitePointTag->data[0].Z;
1280 }
1281 }
1282
1283 /* ViewingConditions: */
1284 {
1285 icmViewingConditions *ro;
1286
1287 if ((ro = (icmViewingConditions *)pp->read_tag(pp, icSigViewingConditionsTag)) != NULL
1288 && ro->ttype == icSigViewingConditionsType) {
1289
1290 /* ro->illuminant.X */
1291 /* ro->illuminant.Z */
1292
1293 Li = ro->illuminant.Y;
1294
1295 /* Reflect illuminant off the media white */
1296 Lvr = Li * Wxyz[1];
1297
1298 /* Illuminant color */
1299 Ixyz[0] = ro->illuminant.X/ro->illuminant.Y;
1300 Ixyz[1] = 1.0;
1301 Ixyz[2] = ro->illuminant.Z/ro->illuminant.Y;
1302
1303 /* Assume ICC surround is CICAM97 background */
1304 /* ro->surround.X */
1305 /* ro->surround.Z */
1306 La = ro->surround.Y;
1307
1308 /* ro->stdIlluminant ie D50, D65, D93, A etc. */
1309 }
1310 }
1311
1312 /* Stuff we might need */
1313
1314 /* Technology: */
1315 {
1316 icmSignature *ro;
1317
1318 /* Try and read the tag from the file */
1319 if ((ro = (icmSignature *)pp->read_tag(pp, icSigTechnologyTag)) != NULL
1320 && ro->ttype != icSigSignatureType) {
1321
1322 tsig = ro->sig;
1323 }
1324 }
1325
1326 devc = pp->header->deviceClass; /* Type of profile */
1327 if (devc == icSigLinkClass
1328 || devc == icSigAbstractClass
1329 || devc == icSigColorSpaceClass
1330 || devc == icSigNamedColorClass)
1331 return 2;
1332
1333 /*
1334 icSigInputClass
1335 icSigDisplayClass
1336 icSigOutputClass
1337 */
1338
1339 if ((pp->header->flags & icTransparency) != 0)
1340 trans = 1;
1341 else
1342 trans = 0;
1343
1344
1345 /* figure Lv if we have the information */
1346 if (Lve >= 0.0)
1347 Lv = Lve; /* Emmisive image white luminance */
1348 else
1349 Lv = Lvr; /* Reflectance image white luminance */
1350
1351 /* Fudge the technology signature */
1352 if (tsig == icMaxEnumTechnology) {
1353 if (devc == icSigDisplayClass)
1354 tsig = icSigCRTDisplay;
1355 }
1356
1357 #ifndef NEVER
1358 printf("Enumeration = %d\n", Ev);
1359 printf("Viewing Conditions:\n");
1360 printf("White adaptation color %f %f %f\n",Wxyz[0], Wxyz[1], Wxyz[2]);
1361 printf("Adapting Luminance La = %f\n",La);
1362 printf("Illuminant color %f %f %f\n",Ixyz[0], Ixyz[1], Ixyz[2]);
1363 printf("Illuminant Luminance Li = %f\n",Li);
1364 printf("Background Luminance Lb = %f\n",Lb);
1365 printf("Relative Background Yb = %f\n",Yb);
1366 printf("Emissive Image White Lve = %f\n",Lve);
1367 printf("Reflective Image White Lvr = %f\n",Lvr);
1368 printf("Device Image White Lv = %f\n",Lv);
1369 printf("Relative Flare Yf = %f\n",Yf);
1370 printf("Relative Glare Yg = %f\n",Yg);
1371 printf("Glare color %f %f %f\n",Gxyz[0], Gxyz[1], Gxyz[2]);
1372 printf("Technology = %s\n",tag2str(tsig));
1373 printf("deviceClass = %s\n",tag2str(devc));
1374 printf("Transparency = %d\n",trans);
1375 // hk ? hkscale ?
1376 #endif
1377
1378 /* See if the viewing conditions are completely defined as ICC can do it */
1379 if (Wxyz[0] >= 0.0 && Wxyz[1] >= 0.0 && Wxyz[2] >= 0.0
1380 && La >= 0.0
1381 && Yb >= 0.0
1382 && Lv >= 0.0
1383 && Yf >= 0.0
1384 && Yg >= 0.0
1385 && Gxyz[0] >= 0.0 && Gxyz[1] >= 0.0 && Gxyz[2] >= 0.0) {
1386
1387 vc->Ev = vc_none;
1388 vc->Wxyz[0] = Wxyz[0];
1389 vc->Wxyz[1] = Wxyz[1];
1390 vc->Wxyz[2] = Wxyz[2];
1391 vc->La = La;
1392 vc->Yb = Yb;
1393 vc->Lv = Lv;
1394 vc->Yf = Yf;
1395 vc->Yg = Yg;
1396 vc->Gxyz[0] = Gxyz[0];
1397 vc->Gxyz[1] = Gxyz[1];
1398 vc->Gxyz[2] = Gxyz[2];
1399 return 0;
1400 }
1401
1402 /* Hmm. We didn't get all the info an ICC can contain. */
1403 /* We will try to guess some reasonable defaults */
1404
1405 /* Have we at least got an adaptation white point ? */
1406 if (Wxyz[0] < 0.0 || Wxyz[1] < 0.0 || Wxyz[2] < 0.0)
1407 return 2; /* No */
1408
1409 /* Have we got the technology ? */
1410 if (tsig == icMaxEnumTechnology)
1411 return 2; /* Hopeless */
1412
1413 /* Guess from the technology */
1414 switch (tsig) {
1415
1416 /* This is inputing either a representation of */
1417 /* a natural scene captured on another a print medium, or */
1418 /* are is assuming that the medium is the original. */
1419 /* We will assume that the print is the original. */
1420 case icSigReflectiveScanner:
1421 {
1422 if (La < 0.0) /* No adapting luminance */
1423 La = 34.0; /* Use a practical print evaluation number */
1424 if (Yb < 0.0) /* No background relative luminance */
1425 Yb = 0.2; /* Assume grey world */
1426 if (Lv < 0.0) /* No device image luminance */
1427 Ev = vc_average; /* Assume average viewing conditions */
1428 if (Yf < 0.0) /* No flare figure */
1429 Yf = 0.0; /* Assume 0% flare */
1430 if (Yg < 0.0) /* No glare figure */
1431 Yg = 0.01; /* Assume 1% glare */
1432 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1433 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1434 break;
1435 }
1436
1437 /* This is inputing either a representation of */
1438 /* a natural scene captured on another a photo medium, or */
1439 /* are is assuming that the medium is the original. */
1440 /* We will assume a compromise media original, natural scene */
1441 case icSigFilmScanner:
1442 {
1443 if (La < 0.0) /* No adapting luminance */
1444 La = 50.0; /* Use bright indoors, dull outdoors */
1445 if (Yb < 0.0) /* No background relative luminance */
1446 Yb = 0.2; /* Assume grey world */
1447 if (Lv < 0.0) /* No device image luminance */
1448 Ev = vc_average; /* Assume average viewing conditions */
1449 if (Yf < 0.0) /* No flare figure */
1450 Yf = 0.0; /* Assume 0% flare */
1451 if (Yg < 0.0) /* No glare figure */
1452 Yg = 0.01; /* Assume 1% glare */
1453 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1454 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1455 break;
1456 }
1457
1458 /* Direct scene to value devices. */
1459 case icSigDigitalCamera:
1460 case icSigVideoCamera:
1461 {
1462 if (La < 0.0) /* No adapting luminance */
1463 La = 110.0; /* Use very bright indoors, usual outdoors */
1464 if (Yb < 0.0) /* No background relative luminance */
1465 Yb = 0.2; /* Assume grey world */
1466 if (Lv < 0.0) /* No device image luminance */
1467 Ev = vc_average; /* Assume average viewing conditions */
1468 if (Yf < 0.0) /* No flare figure */
1469 Yf = 0.0; /* Assume 0% flare */
1470 if (Yg < 0.0) /* No glare figure */
1471 Yg = 0.01; /* Assume 1% glare */
1472 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1473 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1474 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1475 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1476 break;
1477 }
1478
1479 /* Emmisive displays. */
1480 /* Assume a video monitor is in a darker environment than a CRT */
1481 case icSigVideoMonitor:
1482 {
1483 if (La < 0.0) /* No adapting luminance */
1484 La = 4.0; /* Darkened work environment */
1485 if (Yb < 0.0) /* No background relative luminance */
1486 Yb = 0.2; /* Assume grey world */
1487 if (Lv < 0.0) /* No device image luminance */
1488 Ev = vc_dim; /* Assume dim viewing conditions */
1489 if (Yf < 0.0) /* No flare figure */
1490 Yf = 0.0; /* Assume 0% flare */
1491 if (Yg < 0.0) /* No glare figure */
1492 Yg = 0.01; /* Assume 1% glare */
1493 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1494 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1495 break;
1496 }
1497
1498
1499 /* Assume a typical work environment */
1500 case icSigCRTDisplay:
1501 case icSigPMDisplay:
1502 case icSigAMDisplay:
1503 {
1504 if (La < 0.0) /* No adapting luminance */
1505 La = 33.0; /* Typical work environment */
1506 if (Yb < 0.0) /* No background relative luminance */
1507 Yb = 0.2; /* Assume grey world */
1508 if (Lv < 0.0) /* No device image luminance */
1509 Ev = vc_average; /* Assume average viewing conditions */
1510 if (Yf < 0.0) /* No flare figure */
1511 Yf = 0.0; /* Assume 0% flare */
1512 if (Yg < 0.0) /* No glare figure */
1513 Yg = 0.01; /* Assume 1% glare */
1514 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1515 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1516 break;
1517 }
1518
1519 /* Photo CD has its own viewing definitions */
1520 /* (It represents original scene colors) */
1521 case icSigPhotoCD:
1522 {
1523 if (La < 0.0) /* No adapting luminance */
1524 La = 320.0; /* Bright outdoors */
1525 if (Yb < 0.0) /* No background relative luminance */
1526 Yb = 0.2; /* Assume grey world */
1527 if (Lv < 0.0) /* No device image luminance */
1528 Ev = vc_average; /* Assume average viewing conditions */
1529 if (Yf < 0.0) /* No flare figure */
1530 Yf = 0.0; /* Assume 0% flare */
1531 if (Yg < 0.0) /* No glare figure */
1532 Yg = 0.0; /* Assume 0% glare */
1533 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1534 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1535 break;
1536 }
1537
1538 /* Projection devices, either direct, or */
1539 /* via another intermediate medium. */
1540 /* Assume darkened room, little background */
1541 case icSigProjectionTelevision:
1542 {
1543 if (La < 0.0) /* No adapting luminance */
1544 La = 7.0; /* Dark environment */
1545 if (Yb < 0.0) /* No background relative luminance */
1546 Yb = 0.1; /* Assume little background */
1547 if (Lv < 0.0) /* No device image luminance */
1548 Ev = vc_dim; /* Dim environment */
1549 if (Yf < 0.0) /* No flare figure */
1550 Yf = 0.0; /* Assume 0% flare */
1551 if (Yg < 0.0) /* No glare figure */
1552 Yg = 0.01; /* Assume 1% glare */
1553 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1554 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1555 break;
1556 }
1557 /* Assume very darkened room, no background */
1558 case icSigFilmWriter:
1559 {
1560 if (La < 0.0) /* No adapting luminance */
1561 La = 7.0; /* Dark environment */
1562 if (Yb < 0.0) /* No background relative luminance */
1563 Yb = 0.0; /* Assume no background */
1564 if (Lv < 0.0) /* No device image luminance */
1565 Ev = vc_dark; /* Dark environment */
1566 if (Yf < 0.0) /* No flare figure */
1567 Yf = 0.0; /* Assume 0% flare */
1568 if (Yg < 0.0) /* No glare figure */
1569 Yg = 0.01; /* Assume 1% glare */
1570 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1571 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1572 break;
1573 }
1574
1575 /* Printed media devices. */
1576 /* Assume a normal print viewing environment */
1577 case icSigInkJetPrinter:
1578 case icSigThermalWaxPrinter:
1579 case icSigElectrophotographicPrinter:
1580 case icSigElectrostaticPrinter:
1581 case icSigDyeSublimationPrinter:
1582 case icSigPhotographicPaperPrinter:
1583 case icSigPhotoImageSetter:
1584 case icSigGravure:
1585 case icSigOffsetLithography:
1586 case icSigSilkscreen:
1587 case icSigFlexography:
1588 {
1589 if (La < 0.0) /* No adapting luminance */
1590 La = 40.0; /* Use a practical print evaluation number */
1591 if (Yb < 0.0) /* No background relative luminance */
1592 Yb = 0.2; /* Assume grey world */
1593 if (Lv < 0.0) /* No device image luminance */
1594 Ev = vc_average; /* Assume average viewing conditions */
1595 if (Yf < 0.0) /* No flare figure */
1596 Yf = 0.0; /* Assume 0% flare */
1597 if (Yg < 0.0) /* No glare figure */
1598 Yg = 0.01; /* Assume 1% glare */
1599 if (Gxyz[0] < 0.0 || Gxyz[1] < 0.0 || Gxyz[2] < 0.0) /* No flare color */
1600 Gxyz[0] = Wxyz[0], Gxyz[1] = Wxyz[1], Gxyz[2] = Wxyz[2];
1601 break;
1602 }
1603
1604 default:
1605 {
1606 return 2;
1607 }
1608 }
1609
1610 return 1;
1611 }
1612
1613 /* Write our viewing conditions to the underlying ICC profile, */
1614 /* using a private tag. */
xicc_set_viewcond(xicc * p,icxViewCond * vc)1615 void xicc_set_viewcond(
1616 xicc *p, /* Expanded profile we're working with */
1617 icxViewCond *vc /* Viewing parameters to return */
1618 ) {
1619 //icc *pp = p->pp; /* Base ICC */
1620
1621 // ~~1 Not implemented yet
1622 }
1623
1624
1625
1626 /* Return an enumerated viewing condition */
1627 /* Return enumeration if OK, -999 if there is no such enumeration. */
1628 /* xicc may be NULL if just the description is wanted, */
1629 /* or an explicit white point is provided. */
xicc_enum_viewcond(xicc * p,icxViewCond * vc,int no,char * as,int desc,double * wp)1630 int xicc_enum_viewcond(
1631 xicc *p, /* Expanded profile to get white point (May be NULL if desc NZ) */
1632 icxViewCond *vc, /* Viewing parameters to return, May be NULL if desc is nz */
1633 int no, /* Enumeration to return, -1 for default, -2 for none */
1634 char *as, /* String alias to number, NULL if none */
1635 int desc, /* NZ - Just return a description of this enumeration in vc */
1636 double *wp /* Provide XYZ white point if xicc is NULL */
1637 ) {
1638
1639 if (desc == 0) { /* We're setting the viewing condition */
1640 icc *pp; /* Base ICC */
1641 icmXYZArray *whitePointTag;
1642
1643 if (vc == NULL)
1644 return -999;
1645
1646 if (p == NULL) {
1647 if (wp == NULL)
1648 return -999;
1649 vc->Wxyz[0] = wp[0];
1650 vc->Wxyz[1] = wp[1];
1651 vc->Wxyz[2] = wp[2];
1652 } else {
1653
1654 pp = p->pp;
1655 if ((whitePointTag = (icmXYZArray *)pp->read_tag(pp, icSigMediaWhitePointTag)) != NULL
1656 && whitePointTag->ttype == icSigXYZType && whitePointTag->size >= 1) {
1657 vc->Wxyz[0] = whitePointTag->data[0].X;
1658 vc->Wxyz[1] = whitePointTag->data[0].Y;
1659 vc->Wxyz[2] = whitePointTag->data[0].Z;
1660 } else {
1661 if (wp == NULL) {
1662 sprintf(p->err,"Enum VC: Failed to read Media White point");
1663 p->errc = 2;
1664 return -999;
1665 }
1666 vc->Wxyz[0] = wp[0];
1667 vc->Wxyz[1] = wp[1];
1668 vc->Wxyz[2] = wp[2];
1669 }
1670 }
1671
1672 /* Set a default Glare color */
1673 vc->Gxyz[0] = vc->Wxyz[0];
1674 vc->Gxyz[1] = vc->Wxyz[1];
1675 vc->Gxyz[2] = vc->Wxyz[2];
1676
1677 /* Default HK scaling factor = none */
1678 vc->hkscale = 1.0;
1679 }
1680
1681 /*
1682
1683 Typical adapting field luminances and white luminance in reflective setup:
1684 (Note that displays Lv is typically brighter under the same conditions)
1685
1686 E = illuminance in Lux
1687 La = Adapting field luminance in cd/m^2, assuming 20% reflectance from surround
1688 Lv = White luminance assuming 100% reflectance
1689
1690 E La Lv Condition
1691 11 0.7 4 Twilight
1692 32 2 10 Subdued indoor lighting
1693 64 4 20 Less than typical office light; sometimes recommended for
1694 display-only workplaces (sRGB)
1695 350 22 111 Typical Office (sRGB annex D)
1696 500 32 160 Practical print evaluationa (ISO-3664 P2)
1697 1000 64 318 Good Print evaluation (CIE 116-1995)
1698 1000 64 318 Television Studio lighting
1699 1000 64 318 Overcast Outdoors
1700 2000 127 637 Critical print evaluation (ISO-3664 P1)
1701 10000 637 3183 Typical outdoors, full daylight
1702 50000 3185 15915 Bright summers day
1703
1704 Display numbers:
1705
1706 SMPTE video standard white 100
1707 SMPTE cinema standard white 55
1708
1709 Flare is image content dependent, and is typically 1% from factors
1710 including display self illumination and observer/camera internal
1711 stray light. Because image content is not static, using a 1% of white point
1712 flare results quite erronious appearance modelling for predominantly
1713 dark images. As a result, it is best to default to a Yf of 0%,
1714 and only introduce a higher number depending on the known image content.
1715
1716 Glare is assumed to be from the ambient light reflecting from the display
1717 and also striking the observer directly, and is (typically) defaulted
1718 to 1% of ambient here. (too low ? Typical displays are 4-10%)
1719
1720 */
1721
1722 if (no == -1
1723 || (as != NULL && stricmp(as,"d") == 0)) {
1724
1725 no = -1;
1726 if (vc != NULL) {
1727 vc->desc = " d - Default Viewing Condition";
1728 vc->Ev = vc_average; /* Average viewing conditions */
1729 vc->La = 50.0; /* Practical to Good lighting */
1730 vc->Lv = 250.0; /* Average viewing conditions ratio */
1731 vc->Yb = 0.2; /* Grey world */
1732 vc->Yf = 0.0; /* 0% flare */
1733 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1734 }
1735 }
1736 else if (no == 2
1737 || (as != NULL && stricmp(as,"pc") == 0)) {
1738
1739 no = 2;
1740 if (vc != NULL) {
1741 vc->desc = " pc - Critical print evaluation environment (ISO-3664 P1)";
1742 vc->Ev = vc_average; /* Average viewing conditions */
1743 vc->La = 127.0; /* 0.2 * Lv ? */
1744 vc->Lv = 2000.0/3.1415; /* White of the image field */
1745 vc->Yb = 0.2; /* Grey world */
1746 vc->Yf = 0.0; /* 0% flare */
1747 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1748 }
1749 }
1750 else if (no == 0
1751 || (as != NULL && stricmp(as,"pp") == 0)) {
1752
1753 no = 0;
1754 if (vc != NULL) {
1755 vc->desc = " pp - Practical Reflection Print (ISO-3664 P2)";
1756 vc->Ev = vc_none; /* Use explicit La/Lv */
1757 vc->La = 32.0; /* 0.2 * Lv ? */
1758 vc->Lv = 500.0/3.1415; /* White of the image field */
1759 vc->Yb = 0.2; /* Grey world */
1760 vc->Yf = 0.0; /* 0% flare */
1761 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1762 }
1763 }
1764 else if (no == 1
1765 || (as != NULL && stricmp(as,"pe") == 0)) {
1766
1767 no = 1;
1768 if (vc != NULL) {
1769 vc->desc = " pe - Print evaluation environment (CIE 116-1995)";
1770 vc->Ev = vc_none; /* Use explicit La/Lv */
1771 vc->La = 30.0; /* 0.2 * Lv ? */
1772 vc->Lv = 150.0; /* White of the image field */
1773 vc->Yb = 0.2; /* Grey world */
1774 vc->Yf = 0.0; /* 0% flare */
1775 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1776 }
1777 }
1778 else if (no == 4
1779 || (as != NULL && stricmp(as,"mb") == 0)) {
1780
1781 no = 4;
1782 if (vc != NULL) {
1783 vc->desc = " mb - Bright monitor in bright work environment";
1784 vc->Ev = vc_none; /* Use explicit La/Lv */
1785 vc->La = 42.0; /* Bright work environment */
1786 vc->Lv = 150.0; /* White of the image field */
1787 vc->Yb = 0.2; /* Grey world */
1788 vc->Yf = 0.0; /* 0% flare */
1789 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1790 }
1791 }
1792 else if (no == 3
1793 || (as != NULL && stricmp(as,"mt") == 0)) {
1794
1795 no = 3;
1796 if (vc != NULL) {
1797 vc->desc = " mt - Monitor in typical work environment";
1798 vc->Ev = vc_none; /* Use explicit La/Lv */
1799 vc->La = 22.0; /* Typical work environment */
1800 vc->Lv = 120.0; /* White of the image field */
1801 vc->Yb = 0.2; /* Grey world */
1802 vc->Yf = 0.0; /* 0% flare */
1803 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1804 }
1805 }
1806 else if (no == 5
1807 || (as != NULL && stricmp(as,"md") == 0)) {
1808
1809 no = 5;
1810 if (vc != NULL) {
1811 vc->desc = " md - Monitor in darkened work environment";
1812 vc->Ev = vc_none; /* Use explicit La/Lv */
1813 vc->La = 10.0; /* Darkened work environment */
1814 vc->Lv = 100.0; /* White of the image field */
1815 vc->Yb = 0.2; /* Grey world */
1816 vc->Yf = 0.0; /* 0% flare */
1817 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1818 }
1819 }
1820 else if (no == 6
1821 || (as != NULL && stricmp(as,"jm") == 0)) {
1822
1823 no = 6;
1824 if (vc != NULL) {
1825 vc->desc = " jm - Projector in dim environment";
1826 vc->Ev = vc_none; /* Use explicit La/Lv */
1827 vc->La = 10.0; /* Adaptation is from display */
1828 vc->Lv = 80.0; /* White of the image field */
1829 vc->Yb = 0.2; /* Grey world */
1830 vc->Yf = 0.0; /* 0% flare */
1831 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1832 }
1833 }
1834 else if (no == 7
1835 || (as != NULL && stricmp(as,"jd") == 0)) {
1836
1837 no = 7;
1838 if (vc != NULL) {
1839 vc->desc = " jd - Projector in dark environment";
1840 vc->Ev = vc_none; /* Use explicit La/Lv */
1841 vc->La = 8.0; /* Adaptation is from display */
1842 vc->Lv = 80.0; /* White of the image field */
1843 vc->Yb = 0.2; /* Grey world */
1844 vc->Yf = 0.0; /* 0% flare */
1845 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1846 }
1847 }
1848 else if (no == 8
1849 || (as != NULL && stricmp(as,"tv") == 0)) {
1850
1851 no = 8;
1852 if (vc != NULL) {
1853 vc->desc = " tv - Television/Film Studio";
1854 vc->Ev = vc_none; /* Compute from La/Lv */
1855 vc->La = 0.2 * 1000.0/3.1415; /* Adative/Surround */
1856 vc->Yb = 0.2; /* Grey world */
1857 vc->Lv = 1000.0/3.1415; /* White of the image field */
1858 vc->Yf = 0.0; /* 0% flare */
1859 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1860 }
1861 }
1862 else if (no == 9
1863 || (as != NULL && stricmp(as,"pcd") == 0)) {
1864
1865 no = 9;
1866 if (vc != NULL) {
1867 vc->desc = "pcd - Photo CD - original scene outdoors";
1868 vc->Ev = vc_average; /* Average viewing conditions */
1869 vc->La = 320.0; /* Typical outdoors, 1600 cd/m^2 */
1870 vc->Yb = 0.2; /* Grey world */
1871 vc->Yf = 0.0; /* 0% flare */
1872 vc->Yg = 0.0; /* 0% glare - assumed to be compensated ? */
1873 }
1874 }
1875 else if (no == 10
1876 || (as != NULL && stricmp(as,"ob") == 0)) {
1877
1878 no = 10;
1879 if (vc != NULL) {
1880 vc->desc = " ob - Original scene - Bright Outdoors";
1881 vc->Ev = vc_average; /* Average viewing conditions */
1882 vc->La = 2000.0; /* Bright Outdoors */
1883 vc->Yb = 0.2; /* Grey world */
1884 vc->Yf = 0.0; /* 0% flare */
1885 vc->Yg = 0.0; /* 0% glare - assumed to be compensated ? */
1886 }
1887 }
1888 else if (no == 11
1889 || (as != NULL && stricmp(as,"cx") == 0)) {
1890
1891 no = 11;
1892 if (vc != NULL) {
1893 vc->desc = " cx - Cut Sheet Transparencies on a viewing box";
1894 vc->Ev = vc_cut_sheet; /* Cut sheet viewing conditions */
1895 vc->La = 53.0; /* Dim, adapted to slide ? */
1896 vc->Yb = 0.2; /* Grey world */
1897 vc->Yf = 0.0; /* 0% flare */
1898 vc->Yg = 0.01 * XICC_DEFAULT_GLARE; /* 5% glare */
1899 }
1900 }
1901 else {
1902 if (p != NULL) {
1903 sprintf(p->err,"Enum VC: Unrecognised enumeration %d",no);
1904 p->errc = 1;
1905 }
1906 return -999;
1907 }
1908
1909 return no;
1910 }
1911
1912 /* Debug: dump a Viewing Condition to standard out */
xicc_dump_viewcond(icxViewCond * vc)1913 void xicc_dump_viewcond(
1914 icxViewCond *vc
1915 ) {
1916 printf("Viewing Condition:\n");
1917 if (vc->Ev == vc_dark)
1918 printf(" Surround to Image: Dark\n");
1919 else if (vc->Ev == vc_dim)
1920 printf(" Surround to Image: Dim\n");
1921 else if (vc->Ev == vc_average)
1922 printf(" Surround to Image: Average\n");
1923 else if (vc->Ev == vc_cut_sheet)
1924 printf(" Transparency on Light box\n");
1925 printf(" Adapted white = %f %f %f\n",vc->Wxyz[0], vc->Wxyz[1], vc->Wxyz[2]);
1926 printf(" Adapted luminance = %f cd/m^2\n",vc->La);
1927 printf(" Background to image ratio = %f\n",vc->Yb);
1928 if (vc->Ev == vc_none)
1929 printf(" Image luminance = %f cd/m^2\n",vc->Lv);
1930 printf(" Flare to image ratio = %f\n",vc->Yf);
1931 printf(" Glare to adapting/surround ratio = %f\n",vc->Yg);
1932 printf(" Flare color = %f %f %f\n",vc->Gxyz[0], vc->Gxyz[1], vc->Gxyz[2]);
1933 printf(" HK scaling = %f\n",vc->hkscale);
1934 }
1935
1936
1937 /* Debug: dump an Inking setup to standard out */
xicc_dump_inking(icxInk * ik)1938 void xicc_dump_inking(icxInk *ik) {
1939 printf("Inking settings:\n");
1940 if (ik->tlimit < 0.0)
1941 printf("No total limit\n");
1942 else
1943 printf("Total limit = %f%%\n",ik->tlimit * 100.0);
1944
1945 if (ik->klimit < 0.0)
1946 printf("No black limit\n");
1947 else
1948 printf("Black limit = %f%%\n",ik->klimit * 100.0);
1949
1950 if (ik->KonlyLmin)
1951 printf("K only black as locus Lmin\n");
1952 else
1953 printf("Normal black as locus Lmin\n");
1954
1955 if (ik->k_rule == icxKvalue) {
1956 printf("Inking rule is a fixed K target\n");
1957 } if (ik->k_rule == icxKlocus) {
1958 printf("Inking rule is a fixed locus target\n");
1959 } if (ik->k_rule == icxKluma5 || ik->k_rule == icxKluma5k) {
1960 if (ik->k_rule == icxKluma5)
1961 printf("Inking rule is a 5 parameter locus function of L\n");
1962 else
1963 printf("Inking rule is a 5 parameter K function of L\n");
1964 printf("Ksmth = %f\n",ik->c.Ksmth);
1965 printf("Kskew = %f\n",ik->c.Kskew);
1966 printf("Kstle = %f\n",ik->c.Kstle);
1967 printf("Kstpo = %f\n",ik->c.Kstpo);
1968 printf("Kenpo = %f\n",ik->c.Kenpo);
1969 printf("Kenle = %f\n",ik->c.Kenle);
1970 printf("Kshap = %f\n",ik->c.Kshap);
1971 } if (ik->k_rule == icxKl5l || ik->k_rule == icxKl5lk) {
1972 if (ik->k_rule == icxKl5l)
1973 printf("Inking rule is a 2x5 parameter locus function of L and K aux\n");
1974 else
1975 printf("Inking rule is a 2x5 parameter K function of L and K aux\n");
1976 printf("Min Ksmth = %f\n",ik->c.Ksmth);
1977 printf("Min Kskew = %f\n",ik->c.Kskew);
1978 printf("Min Kstle = %f\n",ik->c.Kstle);
1979 printf("Min Kstpo = %f\n",ik->c.Kstpo);
1980 printf("Min Kenpo = %f\n",ik->c.Kenpo);
1981 printf("Min Kenle = %f\n",ik->c.Kenle);
1982 printf("Min Kshap = %f\n",ik->c.Kshap);
1983 printf("Max Ksmth = %f\n",ik->x.Ksmth);
1984 printf("Max Kskew = %f\n",ik->x.Kskew);
1985 printf("Max Kstle = %f\n",ik->x.Kstle);
1986 printf("Max Kstpo = %f\n",ik->x.Kstpo);
1987 printf("Max Kenpo = %f\n",ik->x.Kenpo);
1988 printf("Max Kenle = %f\n",ik->x.Kenle);
1989 printf("Max Kshap = %f\n",ik->x.Kshap);
1990 }
1991 }
1992
1993 /* ------------------------------------------------------ */
1994 /* Gamut Mapping Intent stuff */
1995
1996 /* Return an enumerated gamut mapping intent */
1997 /* Return enumeration if OK, icxIllegalGMIntent if there is no such enumeration. */
xicc_enum_gmapintent(icxGMappingIntent * gmi,int no,char * as)1998 int xicc_enum_gmapintent(
1999 icxGMappingIntent *gmi, /* Gamut Mapping parameters to return */
2000 int no, /* Enumeration selected, icxNoGMIntent for none */
2001 char *as /* Alias string selector, NULL for none */
2002 ) {
2003 #ifdef USE_CAM
2004 int colccas = 0x3; /* Use abs. CAS for abs colorimetric intents */
2005 int perccas = 0x2; /* Use CAS for other intents */
2006 #else
2007 int colccas = 0x1; /* Use abs. Lab for abs colorimetric intents */
2008 int perccas = 0x0; /* Use Lab for other intents */
2009 fprintf(stderr,"!!!!!! Warning, USE_CAM is off in xicc.c !!!!!!\n");
2010 #endif
2011
2012 gmi->hkscale = -1.0; /* Default is to not override viewing condition HK factor */
2013
2014 /* Assert default if no guidance given */
2015 if (no == icxNoGMIntent && as == NULL)
2016 no = icxDefaultGMIntent;
2017
2018 if (no == 0
2019 || no == icxAbsoluteGMIntent
2020 || (as != NULL && stricmp(as,"a") == 0)) {
2021 /* Map Absolute appearance space Jab to Jab and clip out of gamut */
2022 no = 0;
2023 gmi->as = "a";
2024 gmi->desc = " a - Absolute Colorimetric (in Jab) [ICC Absolute Colorimetric]";
2025 gmi->icci = icAbsoluteColorimetric;
2026 gmi->usecas = colccas; /* Use absolute appearance space */
2027 gmi->usemap = 0; /* Don't use gamut mapping */
2028 gmi->greymf = 0.0;
2029 gmi->glumwcpf = 0.0;
2030 gmi->glumwexf = 0.0;
2031 gmi->glumbcpf = 0.0;
2032 gmi->glumbexf = 0.0;
2033 gmi->glumknf = 0.0;
2034 gmi->bph = gmm_noBPadpt; /* No BP adapation */
2035 gmi->gamcpf = 0.0;
2036 gmi->gamexf = 0.0;
2037 gmi->gamcknf = 0.0;
2038 gmi->gamxknf = 0.0;
2039 gmi->gampwf = 0.0;
2040 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2041 gmi->gamswf = 0.0;
2042 gmi->satenh = 0.0; /* No saturation enhancement */
2043 }
2044 else if (no == 1
2045 || (as != NULL && stricmp(as,"aw") == 0)) {
2046
2047 /* I'm not sure how often this intent is useful. It's less likely than */
2048 /* I though that a printer white point won't fit within the gamut */
2049 /* of a display profile, since the display white always has Y = 1.0, */
2050 /* and no paper has better than about 95% reflectance. */
2051 /* Perhaps it may be more useful for targeting printer profiles ? */
2052
2053 /* Map Absolute Jab to Jab and scale source to avoid clipping the white point */
2054 no = 1;
2055 gmi->as = "aw";
2056 gmi->desc = "aw - Absolute Colorimetric (in Jab) with scaling to fit white point";
2057 gmi->icci = icAbsoluteColorimetric;
2058 gmi->usecas = 0x100 | colccas; /* Absolute Appearance space with scaling */
2059 /* to avoid clipping the source white point */
2060 gmi->usemap = 0; /* Don't use gamut mapping */
2061 gmi->greymf = 0.0;
2062 gmi->glumwcpf = 0.0;
2063 gmi->glumwexf = 0.0;
2064 gmi->glumbcpf = 0.0;
2065 gmi->glumbexf = 0.0;
2066 gmi->glumknf = 0.0;
2067 gmi->bph = gmm_noBPadpt; /* No BP adapation */
2068 gmi->gamcpf = 0.0;
2069 gmi->gamexf = 0.0;
2070 gmi->gamcknf = 0.0;
2071 gmi->gamxknf = 0.0;
2072 gmi->gampwf = 0.0;
2073 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2074 gmi->gamswf = 0.0;
2075 gmi->satenh = 0.0; /* No saturation enhancement */
2076 }
2077 else if (no == 2
2078 || (as != NULL && stricmp(as,"aa") == 0)) {
2079
2080 /* Map appearance space Jab to Jab and clip out of gamut */
2081 no = 2;
2082 gmi->as = "aa";
2083 gmi->desc = "aa - Absolute Appearance";
2084 gmi->icci = icRelativeColorimetric;
2085 gmi->usecas = perccas; /* Appearance space */
2086 gmi->usemap = 0; /* Don't use gamut mapping */
2087 gmi->greymf = 0.0;
2088 gmi->glumwcpf = 0.0;
2089 gmi->glumwexf = 0.0;
2090 gmi->glumbcpf = 0.0;
2091 gmi->glumbexf = 0.0;
2092 gmi->glumknf = 0.0;
2093 gmi->bph = gmm_noBPadpt; /* No BP adapation */
2094 gmi->gamcpf = 0.0;
2095 gmi->gamexf = 0.0;
2096 gmi->gamcknf = 0.0;
2097 gmi->gamxknf = 0.0;
2098 gmi->gampwf = 0.0;
2099 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2100 gmi->gamswf = 0.0;
2101 gmi->satenh = 0.0; /* No saturation enhancement */
2102 }
2103 else if (no == 3
2104 || no == icxRelativeGMIntent
2105 || (as != NULL && stricmp(as,"r") == 0)) {
2106
2107 /* Align neutral axes and linearly map white point, then */
2108 /* map appearance space Jab to Jab and clip out of gamut */
2109 no = 3;
2110 gmi->as = "r";
2111 gmi->desc = " r - White Point Matched Appearance [ICC Relative Colorimetric]";
2112 gmi->icci = icRelativeColorimetric;
2113 gmi->usecas = perccas; /* Appearance space */
2114 gmi->usemap = 1; /* Use gamut mapping */
2115 gmi->greymf = 1.0; /* Fully align grey axis */
2116 gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */
2117 gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */
2118 gmi->glumbcpf = 0.0; /* No compression at black end */
2119 gmi->glumbexf = 0.0; /* No expansion at black end */
2120 gmi->glumknf = 0.0;
2121 gmi->bph = gmm_noBPadpt; /* No BP adapation */
2122 gmi->gamcpf = 0.0;
2123 gmi->gamexf = 0.0;
2124 gmi->gamcknf = 0.0;
2125 gmi->gamxknf = 0.0;
2126 gmi->gampwf = 0.0;
2127 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2128 gmi->gamswf = 0.0;
2129 gmi->satenh = 0.0; /* No saturation enhancement */
2130 }
2131 else if (no == 4
2132 || (as != NULL && stricmp(as,"la") == 0)) {
2133
2134 /* Align neutral axes and linearly map white and black points, then */
2135 /* map appearance space Jab to Jab and clip out of gamut */
2136 no = 4;
2137 gmi->as = "la";
2138 gmi->desc = "la - Luminance axis matched Appearance";
2139 gmi->icci = icRelativeColorimetric;
2140 gmi->usecas = perccas; /* Appearance space */
2141 gmi->usemap = 1; /* Use gamut mapping */
2142 gmi->greymf = 1.0; /* Fully align grey axis */
2143 gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */
2144 gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */
2145 gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */
2146 gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */
2147 gmi->glumknf = 0.0; /* No knee on grey mapping */
2148 gmi->bph = gmm_bendBP; /* extent and bend */
2149 gmi->gamcpf = 0.0; /* No gamut compression */
2150 gmi->gamexf = 0.0; /* No gamut expansion */
2151 gmi->gamcknf = 0.0; /* No knee in gamut compress */
2152 gmi->gamxknf = 0.0; /* No knee in gamut expand */
2153 gmi->gampwf = 0.0; /* No Perceptual surface weighting factor */
2154 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2155 gmi->gamswf = 0.0; /* No Saturation surface weighting factor */
2156 gmi->satenh = 0.0; /* No saturation enhancement */
2157 }
2158 else if (no == 5
2159 || no == icxDefaultGMIntent
2160 || no == icxPerceptualGMIntent
2161 || (as != NULL && stricmp(as,"p") == 0)) {
2162
2163 /* Align neutral axes and perceptually map white and black points, */
2164 /* perceptually compress out of gamut and map appearance space Jab to Jab. */
2165 no = 5;
2166 gmi->as = "p";
2167 gmi->desc = " p - Perceptual (Preferred) (Default) [ICC Perceptual]";
2168 gmi->icci = icPerceptual;
2169 gmi->usecas = perccas; /* Appearance space */
2170 gmi->usemap = 1; /* Use gamut mapping */
2171 gmi->greymf = 1.0; /* Fully align grey axis */
2172 gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */
2173 gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */
2174 gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */
2175 gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */
2176 gmi->glumknf = 1.0; /* Sigma knee in grey compress/expand */
2177 gmi->bph = gmm_bendBP; /* extent and bend */
2178 gmi->gamcpf = 1.0; /* Full gamut compression */
2179 gmi->gamexf = 0.0; /* No gamut expansion */
2180 gmi->gamcknf = 1.0; /* Full Sigma knee in gamut compress */
2181 gmi->gamxknf = 0.0; /* No knee in gamut expand */
2182 gmi->gampwf = 1.0; /* Full Perceptual surface weighting factor */
2183 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2184 gmi->gamswf = 0.0; /* No Saturation surface weighting factor */
2185 gmi->satenh = 0.0; /* No saturation enhancement */
2186 }
2187 else if (no == 6
2188 || (as != NULL && stricmp(as,"pa") == 0)) {
2189
2190 /* Don't align neutral axes, but perceptually compress out of gamut */
2191 /* and map appearance space Jab to Jab. */
2192 no = 6;
2193 gmi->as = "pa";
2194 gmi->desc = "pa - Perceptual Apperance ";
2195 gmi->icci = icPerceptual;
2196 gmi->usecas = perccas; /* Appearance space */
2197 gmi->usemap = 1; /* Use gamut mapping */
2198 gmi->greymf = 0.0; /* Don't align grey axis */
2199 gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */
2200 gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */
2201 gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */
2202 gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */
2203 gmi->glumknf = 1.0; /* Sigma knee in grey compress/expand */
2204 gmi->bph = gmm_bendBP; /* extent and bend */
2205 gmi->gamcpf = 1.0; /* Full gamut compression */
2206 gmi->gamexf = 0.0; /* No gamut expansion */
2207 gmi->gamcknf = 1.0; /* Full Sigma knee in gamut compress */
2208 gmi->gamxknf = 0.0; /* No knee in gamut expand */
2209 gmi->gampwf = 1.0; /* Full Perceptual surface weighting factor */
2210 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2211 gmi->gamswf = 0.0; /* No Saturation surface weighting factor */
2212 gmi->satenh = 0.0; /* No saturation enhancement */
2213 }
2214 else if (no == 7
2215 || (as != NULL && stricmp(as,"lp") == 0)) {
2216
2217 /* Align neutral axes and perceptually map white and black points, */
2218 /* perceptually compress out of gamut and map appearance space Jab to Jab, */
2219 /* and heavily weight preserving the luminance over saturation. */
2220 /* No neutral axis sigma enhancement. */
2221 no = 7;
2222 gmi->as = "lp";
2223 gmi->desc = "lp - Luminance Preserving Perceptual";
2224 gmi->icci = icPerceptual;
2225 gmi->usecas = perccas; /* Appearance space */
2226 gmi->usemap = 1; /* Use gamut mapping */
2227 gmi->greymf = 1.0; /* Fully align grey axis */
2228 gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */
2229 gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */
2230 gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */
2231 gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */
2232 gmi->glumknf = 0.3; /* Low Sigma knee in grey compress/expand */
2233 gmi->bph = gmm_bendBP; /* extent and bend */
2234 gmi->gamcpf = 1.0; /* Full gamut compression */
2235 gmi->gamexf = 0.0; /* No gamut expansion */
2236 gmi->gamcknf = 1.3; /* [1.3] High Sigma knee in gamut compress */
2237 gmi->gamxknf = 0.0; /* No knee in gamut expand */
2238 gmi->gampwf = 0.0; /* No Perceptual weighting factor */
2239 gmi->gamlpwf = 1.0; /* Full Linear Preserving Perceptual wghtg. factor */
2240 gmi->gamswf = 0.0; /* No Saturation weighting factor */
2241 gmi->satenh = 0.0; /* No saturation enhancement */
2242 gmi->hkscale = 0.2; /* Mostly disable HK appearance modeling */
2243 }
2244 else if (no == 8
2245 || (as != NULL && stricmp(as,"ms") == 0)) {
2246
2247 /* Align neutral axes and perceptually map white and black points, */
2248 /* perceptually compress and expand to match gamuts and map Jab to Jab. */
2249 no = 8;
2250 gmi->as = "ms";
2251 gmi->desc = "ms - Saturation";
2252 gmi->icci = icSaturation;
2253 gmi->usecas = perccas; /* Appearance space */
2254 gmi->usemap = 1; /* Use gamut mapping */
2255 gmi->greymf = 1.0; /* Fully align grey axis */
2256 gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */
2257 gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */
2258 gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */
2259 gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */
2260 gmi->glumknf = 1.0; /* Sigma knee in grey compress/expand */
2261 gmi->bph = gmm_bendBP; /* extent and bend */
2262 gmi->gamcpf = 1.0; /* Full gamut compression */
2263 gmi->gamexf = 1.0; /* Full gamut expansion */
2264 gmi->gamcknf = 1.1; /* Sigma knee in gamut compress */
2265 gmi->gamxknf = 0.4; /* Moderate Sigma knee in gamut expand */
2266 gmi->gampwf = 0.2; /* Slight perceptual surface weighting factor */
2267 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2268 gmi->gamswf = 0.8; /* Most saturation surface weighting factor */
2269 gmi->satenh = 0.0; /* No saturation enhancement */
2270 }
2271 else if (no == 9
2272 || no == icxSaturationGMIntent
2273 || (as != NULL && stricmp(as,"s") == 0)) {
2274
2275 /* Same as "ms" but enhance saturation */
2276 no = 9;
2277 gmi->as = "s";
2278 gmi->desc = " s - Enhanced Saturation [ICC Saturation]";
2279 gmi->icci = icSaturation;
2280 gmi->usecas = perccas; /* Appearance space */
2281 gmi->usemap = 1; /* Use gamut mapping */
2282 gmi->greymf = 1.0; /* Fully align grey axis */
2283 gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */
2284 gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */
2285 gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */
2286 gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */
2287 gmi->glumknf = 1.0; /* Sigma knee in grey compress/expand */
2288 gmi->bph = gmm_bendBP; /* extent and bend */
2289 gmi->gamcpf = 1.0; /* Full gamut compression */
2290 gmi->gamexf = 1.0; /* Full gamut expansion */
2291 gmi->gamcknf = 1.1; /* High sigma knee in gamut compress */
2292 gmi->gamxknf = 0.5; /* Moderate sigma knee in gamut expand */
2293 gmi->gampwf = 0.0; /* No Perceptual surface weighting factor */
2294 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2295 gmi->gamswf = 1.0; /* Full Saturation surface weighting factor */
2296 gmi->satenh = 0.9; /* Medium saturation enhancement */
2297 }
2298 else if (no == 10
2299 || (as != NULL && stricmp(as,"al") == 0)) {
2300
2301 /* Map absolute L*a*b* to L*a*b* and clip out of gamut */
2302 no = 10;
2303 gmi->as = "al";
2304 gmi->desc = "al - Absolute Colorimetric (Lab)";
2305 gmi->icci = icAbsoluteColorimetric;
2306 gmi->usecas = 0x1; /* Don't use appearance space, use abs. L*a*b* */
2307 gmi->usemap = 0; /* Don't use gamut mapping */
2308 gmi->greymf = 0.0;
2309 gmi->glumwcpf = 0.0;
2310 gmi->glumwexf = 0.0;
2311 gmi->glumbcpf = 0.0;
2312 gmi->glumbexf = 0.0;
2313 gmi->glumknf = 0.0;
2314 gmi->bph = gmm_noBPadpt; /* No BP adapation */
2315 gmi->gamcpf = 0.0;
2316 gmi->gamexf = 0.0;
2317 gmi->gamcknf = 0.0;
2318 gmi->gamxknf = 0.0;
2319 gmi->gampwf = 0.0;
2320 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2321 gmi->gamswf = 0.0;
2322 gmi->satenh = 0.0; /* No saturation enhancement */
2323 }
2324 else if (no == 11
2325 || (as != NULL && stricmp(as,"rl") == 0)) {
2326
2327 /* Align neutral axes and linearly map white point, then */
2328 /* map L*a*b* to L*a*b* and clip out of gamut */
2329 no = 11;
2330 gmi->as = "rl";
2331 gmi->desc = "rl - White Point Matched Colorimetric (Lab)";
2332 gmi->icci = icRelativeColorimetric;
2333 gmi->usecas = 0x0; /* Don't use appearance space, use relative L*a*b* */
2334 gmi->usemap = 1; /* Use gamut mapping */
2335 gmi->greymf = 1.0; /* And linearly map white point */
2336 gmi->glumwcpf = 1.0;
2337 gmi->glumwexf = 1.0;
2338 gmi->glumbcpf = 0.0;
2339 gmi->glumbexf = 0.0;
2340 gmi->glumknf = 0.0;
2341 gmi->bph = gmm_noBPadpt; /* No BP adapation */
2342 gmi->gamcpf = 0.0;
2343 gmi->gamexf = 0.0;
2344 gmi->gamcknf = 0.0;
2345 gmi->gamxknf = 0.0;
2346 gmi->gampwf = 0.0;
2347 gmi->gamlpwf = 0.0; /* No Linear Preserving Perceptual surface wghtg. factor */
2348 gmi->gamswf = 0.0;
2349 gmi->satenh = 0.0; /* No saturation enhancement */
2350 }
2351 else { /* icxIllegalGMIntent */
2352 return icxIllegalGMIntent;
2353 }
2354
2355 return no;
2356 }
2357
2358
2359 /* Debug: dump a Gamut Mapping specification */
xicc_dump_gmi(icxGMappingIntent * gmi)2360 void xicc_dump_gmi(
2361 icxGMappingIntent *gmi /* Gamut Mapping parameters to return */
2362 ) {
2363 printf(" Gamut Mapping Specification:\n");
2364 if (gmi->desc != NULL)
2365 printf(" Description = '%s'\n",gmi->desc);
2366 printf(" Closest ICC intent = '%s'\n",icm2str(icmRenderingIntent,gmi->icci));
2367
2368 if ((gmi->usecas & 0xff) == 0)
2369 printf(" Not using Color Apperance Space - using L*a*b*\n");
2370 else if ((gmi->usecas & 0xff) == 1)
2371 printf(" Not using Color Apperance Space - using Absoute L*a*b*\n");
2372 else if ((gmi->usecas & 0xff) == 2)
2373 printf(" Using Color Apperance Space\n");
2374 else if ((gmi->usecas & 0xff) == 3)
2375 printf(" Using Absolute Color Apperance Space\n");
2376
2377 if ((gmi->usecas & 0x100) != 0)
2378 printf(" Scaling source to avoid white point clipping\n");
2379
2380 if (gmi->usemap == 0)
2381 printf(" Not using Mapping\n");
2382 else {
2383 printf(" Using Mapping with parameters:\n");
2384 printf(" Grey axis alignment factor %f\n", gmi->greymf);
2385 printf(" Grey axis white compression factor %f\n", gmi->glumwcpf);
2386 printf(" Grey axis white expansion factor %f\n", gmi->glumwexf);
2387 printf(" Grey axis black compression factor %f\n", gmi->glumbcpf);
2388 printf(" Grey axis black expansion factor %f\n", gmi->glumbexf);
2389 printf(" Grey axis knee factor %f\n", gmi->glumknf);
2390 printf(" Black point algorithm: ");
2391 switch(gmi->bph) {
2392 case gmm_clipBP: printf("Neutral axis no-adapt extend and clip\n"); break;
2393 case gmm_BPadpt: printf("Neutral axis fully adapt\n"); break;
2394 case gmm_bendBP: printf("Neutral axis no-adapt extend and bend\n"); break;
2395 case gmm_noBPadpt: printf("Neutral axis no-adapt\n"); break;
2396 }
2397 printf(" Gamut compression factor %f\n", gmi->gamcpf);
2398 printf(" Gamut expansion factor %f\n", gmi->gamexf);
2399 printf(" Gamut compression knee factor %f\n", gmi->gamcknf);
2400 printf(" Gamut expansion knee factor %f\n", gmi->gamxknf);
2401 printf(" Gamut Perceptual mapping weighting factor %f\n", gmi->gampwf);
2402 printf(" Gamut Lightness Preserving Perceptual mapping weighting %f\n", gmi->gamlpwf);
2403 printf(" Gamut Saturation mapping weighting factor %f\n", gmi->gamswf);
2404 printf(" Saturation enhancement factor %f\n", gmi->satenh);
2405 }
2406 if (gmi->hkscale >= 0.0)
2407 printf(" HK scale override %f\n", gmi->hkscale);
2408 }
2409
2410 /* ------------------------------------------------------ */
2411 /* Turn xicc xcal into limit calibration callback */
2412
2413 /* Given an icc profile, try and create an xcal */
2414 /* Return NULL on error or no cal */
xiccReadCalTag(icc * p)2415 xcal *xiccReadCalTag(icc *p) {
2416 xcal *cal = NULL;
2417 icTagSignature sig = icmMakeTag('t','a','r','g');
2418 icmText *ro;
2419 int oi, tab;
2420
2421 //printf("~1 about to look for CAL in profile\n");
2422 if ((ro = (icmText *)p->read_tag(p, sig)) != NULL) {
2423 cgatsFile *cgf;
2424 cgats *icg;
2425
2426 if (ro->ttype != icSigTextType)
2427 return NULL;
2428
2429 //printf("~1 found 'targ' tag\n");
2430 if ((icg = new_cgats()) == NULL) {
2431 return NULL;
2432 }
2433 if ((cgf = new_cgatsFileMem(ro->data, ro->size)) != NULL) {
2434 icg->add_other(icg, "CTI3");
2435 oi = icg->add_other(icg, "CAL");
2436
2437 //printf("~1 created cgats object from 'targ' tag\n");
2438 if (icg->read(icg, cgf) == 0) {
2439
2440 for (tab = 0; tab < icg->ntables; tab++) {
2441 if (icg->t[tab].tt == tt_other && icg->t[tab].oi == oi) {
2442 break;
2443 }
2444 }
2445 if (tab < icg->ntables) {
2446 //printf("~1 found CAL table\n");
2447
2448 if ((cal = new_xcal()) == NULL) {
2449 icg->del(icg);
2450 cgf->del(cgf);
2451 return NULL;
2452 }
2453 if (cal->read_cgats(cal, icg, tab, "'targ' tag") != 0) {
2454 #ifdef DEBUG
2455 printf("read_cgats on cal tag failed\n");
2456 #endif
2457 cal->del(cal);
2458 cal = NULL;
2459 }
2460 //else printf("~1 read CAL and creaded xcal object OK\n");
2461 }
2462 }
2463 cgf->del(cgf);
2464 }
2465 icg->del(icg);
2466 }
2467 return cal;
2468 }
2469
2470 /* A callback that uses an xcal, that can be used with icc get_tac */
xiccCalCallback(void * cntx,double * out,double * in)2471 void xiccCalCallback(void *cntx, double *out, double *in) {
2472 xcal *cal = (xcal *)cntx;
2473
2474 cal->interp(cal, out, in);
2475 }
2476
2477 /* ---------------------------------------------- */
2478
2479 /* Utility function - given an open icc profile, */
2480 /* guess which channel is the black. */
2481 /* Return -1 if there is no black channel or it can't be guessed */
icxGuessBlackChan(icc * p)2482 int icxGuessBlackChan(icc *p) {
2483 int kch = -1;
2484
2485 switch (p->header->colorSpace) {
2486 case icSigCmykData:
2487 kch = 3;
2488 break;
2489
2490 /* Use a heuristic to detect the black channel. */
2491 /* This duplicates code in icxLu_comp_bk_point() :-( */
2492 /* Colorant guessing should go in icclib ? */
2493 case icSig2colorData:
2494 case icSig3colorData:
2495 case icSig4colorData:
2496 case icSig5colorData:
2497 case icSig6colorData:
2498 case icSig7colorData:
2499 case icSig8colorData:
2500 case icSig9colorData:
2501 case icSig10colorData:
2502 case icSig11colorData:
2503 case icSig12colorData:
2504 case icSig13colorData:
2505 case icSig14colorData:
2506 case icSig15colorData:
2507 case icSigMch5Data:
2508 case icSigMch6Data:
2509 case icSigMch7Data:
2510 case icSigMch8Data: {
2511 icmLuBase *lu;
2512 double dval[MAX_CHAN];
2513 double ncval[3];
2514 double cvals[MAX_CHAN][3];
2515 int inn, e, nlighter, ndarker;
2516
2517 /* Grab a lookup object */
2518 if ((lu = p->get_luobj(p, icmFwd, icRelativeColorimetric, icSigLabData, icmLuOrdNorm)) == NULL)
2519 error("icxGetLimits: assert: getting Fwd Lookup failed!");
2520
2521 lu->spaces(lu, NULL, &inn, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
2522
2523 /* Decide if the colorspace is aditive or subtractive */
2524
2525 /* First the no colorant value */
2526 for (e = 0; e < inn; e++)
2527 dval[e] = 0.0;
2528 lu->lookup(lu, ncval, dval);
2529
2530 /* Then all the colorants */
2531 nlighter = ndarker = 0;
2532 for (e = 0; e < inn; e++) {
2533 dval[e] = 1.0;
2534 lu->lookup(lu, cvals[e], dval);
2535 dval[e] = 0.0;
2536 if (fabs(cvals[e][0] - ncval[0]) > 5.0) {
2537 if (cvals[e][0] > ncval[0])
2538 nlighter++;
2539 else
2540 ndarker++;
2541 }
2542 }
2543
2544 if (ndarker > 0 && nlighter == 0) { /* Assume subtractive. */
2545 double pbk[3] = { 0.0,0.0,0.0 }; /* Perfect black */
2546 double smd = 1e10; /* Smallest distance */
2547
2548 /* Guess the black channel */
2549 for (e = 0; e < inn; e++) {
2550 double tt;
2551 tt = icmNorm33sq(pbk, cvals[e]);
2552 if (tt < smd) {
2553 smd = tt;
2554 kch = e;
2555 }
2556 }
2557 /* See if the black seems sane */
2558 if (cvals[kch][0] > 40.0
2559 || fabs(cvals[kch][1]) > 10.0
2560 || fabs(cvals[kch][2]) > 10.0) {
2561 kch = -1;
2562 }
2563 }
2564 lu->del(lu);
2565 }
2566 break;
2567
2568 default:
2569 break;
2570 }
2571
2572 return kch;
2573 }
2574
2575 /* Utility function - given an open icc profile, */
2576 /* estmate the total ink limit and black ink limit. */
2577 /* Note that this is rather rough, because ICC profiles */
2578 /* don't have a tag for this information, and ICC profiles */
2579 /* don't have any straightforward way of identifying particular */
2580 /* color channels for > 4 color. */
2581 /* If there are no limits, or they are not discoverable or */
2582 /* applicable, return values of -1.0 */
2583
icxGetLimits(xicc * xp,double * tlimit,double * klimit)2584 void icxGetLimits(
2585 xicc *xp,
2586 double *tlimit,
2587 double *klimit
2588 ) {
2589 icc *p = xp->pp;
2590 int nch;
2591 double max[MAX_CHAN]; /* Max of each channel */
2592 double total;
2593
2594 total = p->get_tac(p, max, xp->cal != NULL ? xiccCalCallback : NULL, (void *)xp->cal);
2595
2596 if (total < 0.0) { /* Not valid */
2597 if (tlimit != NULL)
2598 *tlimit = -1.0;
2599 if (klimit != NULL)
2600 *klimit = -1.0;
2601 return;
2602 }
2603
2604 nch = icmCSSig2nchan(p->header->colorSpace);
2605
2606 /* No effective limit */
2607 if (tlimit != NULL) {
2608 if (total >= (double)nch) {
2609 *tlimit = -1.0;
2610 } else {
2611 *tlimit = total;
2612 }
2613 }
2614
2615 if (klimit != NULL) {
2616 int kch;
2617
2618 kch = icxGuessBlackChan(p);
2619
2620 if (kch < 0 || max[kch] >= 1.0) {
2621 *klimit = -1.0;
2622 } else {
2623 *klimit = max[kch];
2624 }
2625 }
2626 }
2627
2628 /* Replace a non-set limit (ie. < 0.0) with the heuristic from */
2629 /* the given profile. */
icxDefaultLimits(xicc * xp,double * tlout,double tlin,double * klout,double klin)2630 void icxDefaultLimits(
2631 xicc *xp,
2632 double *tlout,
2633 double tlin,
2634 double *klout,
2635 double klin
2636 ) {
2637 if (tlin < 0.0 || klin < 0.0) {
2638 double tl, kl;
2639
2640 icxGetLimits(xp, &tl, &kl);
2641
2642 if (tlin < 0.0)
2643 tlin = tl;
2644
2645 if (klin < 0.0)
2646 klin = kl;
2647 }
2648
2649 if (tlout != NULL)
2650 *tlout = tlin;
2651
2652 if (klout != NULL)
2653 *klout = klin;
2654 }
2655
2656 /* Structure to hold optimisation information */
2657 typedef struct {
2658 xcal *cal;
2659 double ilimit;
2660 double uilimit;
2661 } ulimctx;
2662
2663 /* Callback to find equivalent underlying total limit */
2664 /* and try and maximize it while remaining within gamut */
ulimitfunc(void * cntx,double pv[])2665 static double ulimitfunc(void *cntx, double pv[]) {
2666 ulimctx *cx = (ulimctx *)cntx;
2667 xcal *cal = cx->cal;
2668 int devchan = cal->devchan;
2669 int i;
2670 double dv, odv;
2671 double og = 0.0, rv = 0.0;
2672
2673 double usum = 0.0, sum = 0.0;
2674
2675 /* Comute calibrated sum of channels except last */
2676 for (i = 0; i < (devchan-1); i++) {
2677 double dv = pv[i]; /* Underlying (pre-calibration) device value */
2678 usum += dv; /* Underlying sum */
2679 if (dv < 0.0) {
2680 og += -dv;
2681 dv = 0.0;
2682 } else if (dv > 1.0) {
2683 og += dv - 1.0;
2684 dv = 1.0;
2685 } else
2686 dv = cal->interp_ch(cal, i, dv); /* Calibrated device value */
2687 sum += dv; /* Calibrated device sum */
2688 }
2689 /* Compute the omitted channel value */
2690 dv = cx->ilimit - sum; /* Omitted calibrated device value */
2691 if (dv < 0.0) {
2692 og += -dv;
2693 dv = 0.0;
2694 } else if (dv > 1.0) {
2695 og += dv - 1.0;
2696 dv = 1.0;
2697 } else
2698 dv = cal->inv_interp_ch(cal, i, dv); /* Omitted underlying device value */
2699 usum += dv; /* Underlying sum */
2700 cx->uilimit = usum;
2701
2702 rv = 10000.0 * og - usum; /* Penalize out of gamut, maximize underlying sum */
2703
2704 //printf("~1 returning %f from %f %f %f %f\n",rv,pv[0],pv[1],pv[2],dv);
2705 return rv;
2706 }
2707
2708 /* Given a calibrated total ink limit and an xcal, return the */
2709 /* equivalent underlying (pre-calibration) total ink limit. */
2710 /* This is the maximum equivalent, that makes sure that */
2711 /* the calibrated limit is met or exceeded. */
icxMaxUnderlyingLimit(xcal * cal,double ilimit)2712 double icxMaxUnderlyingLimit(xcal *cal, double ilimit) {
2713 ulimctx cx;
2714 int i;
2715 double dv[MAX_CHAN];
2716 double sr[MAX_CHAN];
2717 double rv; /* Residual value */
2718
2719 if (cal->devchan <= 1) {
2720 return cal->inv_interp_ch(cal, 0, ilimit);
2721 }
2722
2723 cx.cal = cal;
2724 cx.ilimit = ilimit;
2725
2726 for (i = 0; i < (cal->devchan-1); i++) {
2727 sr[i] = 0.05;
2728 dv[i] = 0.1;
2729 }
2730 if (powell(&rv, cal->devchan-1, dv, sr, 0.000001, 1000, ulimitfunc,
2731 (void *)&cx, NULL, NULL) != 0) {
2732 warning("icxUnderlyingLimit() failed for chan %d, ilimit %f\n",cal->devchan,ilimit);
2733 return ilimit;
2734 }
2735 ulimitfunc((void *)&cx, dv);
2736
2737 return cx.uilimit;
2738 }
2739
2740 /* ------------------------------------------------------ */
2741 /* Conversion and deltaE formular that include partial */
2742 /* derivatives, for use within fit parameter optimisations. */
2743
2744 /* CIE XYZ to perceptual Lab with partial derivatives. */
icxdXYZ2Lab(icmXYZNumber * w,double * out,double dout[3][3],double * in)2745 void icxdXYZ2Lab(icmXYZNumber *w, double *out, double dout[3][3], double *in) {
2746 double wp[3], tin[3], dtin[3];
2747 int i;
2748
2749 wp[0] = w->X, wp[1] = w->Y, wp[2] = w->Z;
2750
2751 for (i = 0; i < 3; i++) {
2752 tin[i] = in[i]/wp[i];
2753 dtin[i] = 1.0/wp[i];
2754
2755 if (tin[i] > 0.008856451586) {
2756 dtin[i] *= pow(tin[i], -2.0/3.0) / 3.0;
2757 tin[i] = pow(tin[i],1.0/3.0);
2758 } else {
2759 dtin[i] *= 7.787036979;
2760 tin[i] = 7.787036979 * tin[i] + 16.0/116.0;
2761 }
2762 }
2763
2764 out[0] = 116.0 * tin[1] - 16.0;
2765 dout[0][0] = 0.0;
2766 dout[0][1] = 116.0 * dtin[1];
2767 dout[0][2] = 0.0;
2768
2769 out[1] = 500.0 * (tin[0] - tin[1]);
2770 dout[1][0] = 500.0 * dtin[0];
2771 dout[1][1] = 500.0 * -dtin[1];
2772 dout[1][2] = 0.0;
2773
2774 out[2] = 200.0 * (tin[1] - tin[2]);
2775 dout[2][0] = 0.0;
2776 dout[2][1] = 200.0 * dtin[1];
2777 dout[2][2] = 200.0 * -dtin[2];
2778 }
2779
2780
2781 /* Return the normal Delta E squared, given two Lab values, */
2782 /* including partial derivatives. */
icxdLabDEsq(double dout[2][3],double * Lab0,double * Lab1)2783 double icxdLabDEsq(double dout[2][3], double *Lab0, double *Lab1) {
2784 double rv = 0.0, tt;
2785
2786 tt = Lab0[0] - Lab1[0];
2787 dout[0][0] = 2.0 * tt;
2788 dout[1][0] = -2.0 * tt;
2789 rv += tt * tt;
2790 tt = Lab0[1] - Lab1[1];
2791 dout[0][1] = 2.0 * tt;
2792 dout[1][1] = -2.0 * tt;
2793 rv += tt * tt;
2794 tt = Lab0[2] - Lab1[2];
2795 dout[0][2] = 2.0 * tt;
2796 dout[1][2] = -2.0 * tt;
2797 rv += tt * tt;
2798 return rv;
2799 }
2800
2801 /* Return the CIE94 Delta E color difference measure, squared */
2802 /* including partial derivatives. */
icxdCIE94sq(double dout[2][3],double Lab0[3],double Lab1[3])2803 double icxdCIE94sq(double dout[2][3], double Lab0[3], double Lab1[3]) {
2804 double desq, _desq[2][3];
2805 double dlsq;
2806 double dcsq, _dcsq[2][2]; /* == [x][1,2] */
2807 double c12, _c12[2][2]; /* == [x][1,2] */
2808 double dhsq, _dhsq[2][2]; /* == [x][1,2] */
2809 double rv;
2810
2811 {
2812 double dl, da, db;
2813
2814 dl = Lab0[0] - Lab1[0];
2815 dlsq = dl * dl; /* dl squared */
2816 da = Lab0[1] - Lab1[1];
2817 db = Lab0[2] - Lab1[2];
2818
2819
2820 /* Compute normal Lab delta E squared */
2821 desq = dlsq + da * da + db * db;
2822 _desq[0][0] = 2.0 * dl;
2823 _desq[1][0] = -2.0 * dl;
2824 _desq[0][1] = 2.0 * da;
2825 _desq[1][1] = -2.0 * da;
2826 _desq[0][2] = 2.0 * db;
2827 _desq[1][2] = -2.0 * db;
2828 }
2829
2830 {
2831 double c1, c2, dc, tt;
2832
2833 /* Compute chromanance for the two colors */
2834 c1 = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]);
2835 c2 = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]);
2836 c12 = sqrt(c1 * c2); /* Symetric chromanance */
2837
2838 tt = 0.5 * (pow(c2, 0.5) + 1e-12)/(pow(c1, 1.5) + 1e-12);
2839 _c12[0][0] = Lab0[1] * tt;
2840 _c12[0][1] = Lab0[2] * tt;
2841 tt = 0.5 * (pow(c1, 0.5) + 1e-12)/(pow(c2, 1.5) + 1e-12);
2842 _c12[1][0] = Lab1[1] * tt;
2843 _c12[1][1] = Lab1[2] * tt;
2844
2845 /* delta chromanance squared */
2846 dc = c2 - c1;
2847 dcsq = dc * dc;
2848 if (c1 < 1e-12 || c2 < 1e-12) {
2849 c1 += 1e-12;
2850 c2 += 1e-12;
2851 }
2852 _dcsq[0][0] = -2.0 * Lab0[1] * (c2 - c1)/c1;
2853 _dcsq[0][1] = -2.0 * Lab0[2] * (c2 - c1)/c1;
2854 _dcsq[1][0] = 2.0 * Lab1[1] * (c2 - c1)/c2;
2855 _dcsq[1][1] = 2.0 * Lab1[2] * (c2 - c1)/c2;
2856 }
2857
2858 /* Compute delta hue squared */
2859 dhsq = desq - dlsq - dcsq;
2860 if (dhsq >= 0.0) {
2861 _dhsq[0][0] = _desq[0][1] - _dcsq[0][0];
2862 _dhsq[0][1] = _desq[0][2] - _dcsq[0][1];
2863 _dhsq[1][0] = _desq[1][1] - _dcsq[1][0];
2864 _dhsq[1][1] = _desq[1][2] - _dcsq[1][1];
2865 } else {
2866 dhsq = 0.0;
2867 _dhsq[0][0] = 0.0;
2868 _dhsq[0][1] = 0.0;
2869 _dhsq[1][0] = 0.0;
2870 _dhsq[1][1] = 0.0;
2871 }
2872
2873 {
2874 double sc, scsq, scf;
2875 double sh, shsq, shf;
2876
2877 /* Weighting factors for delta chromanance & delta hue */
2878 sc = 1.0 + 0.048 * c12;
2879 scsq = sc * sc;
2880
2881 sh = 1.0 + 0.014 * c12;
2882 shsq = sh * sh;
2883
2884 rv = dlsq + dcsq/scsq + dhsq/shsq;
2885
2886 scf = 0.048 * -2.0 * dcsq/(scsq * sc);
2887 shf = 0.014 * -2.0 * dhsq/(shsq * sh);
2888 dout[0][0] = _desq[0][0];
2889 dout[0][1] = _dcsq[0][0]/scsq + _c12[0][0] * scf
2890 + _dhsq[0][0]/shsq + _c12[0][0] * shf;
2891 dout[0][2] = _dcsq[0][1]/scsq + _c12[0][1] * scf
2892 + _dhsq[0][1]/shsq + _c12[0][1] * shf;
2893 dout[1][0] = _desq[1][0];
2894 dout[1][1] = _dcsq[1][0]/scsq + _c12[1][0] * scf
2895 + _dhsq[1][0]/shsq + _c12[1][0] * shf;
2896 dout[1][2] = _dcsq[1][1]/scsq + _c12[1][1] * scf
2897 + _dhsq[1][1]/shsq + _c12[1][1] * shf;
2898 return rv;
2899 }
2900 }
2901
2902 // ~~99 not sure if these are correct:
2903
2904 /* Return the normal Delta E given two Lab values, */
2905 /* including partial derivatives. */
icxdLabDE(double dout[2][3],double * Lab0,double * Lab1)2906 double icxdLabDE(double dout[2][3], double *Lab0, double *Lab1) {
2907 double rv = 0.0, tt;
2908
2909 tt = Lab0[0] - Lab1[0];
2910 dout[0][0] = 1.0 * tt;
2911 dout[1][0] = -1.0 * tt;
2912 rv += tt * tt;
2913 tt = Lab0[1] - Lab1[1];
2914 dout[0][1] = 1.0 * tt;
2915 dout[1][1] = -1.0 * tt;
2916 rv += tt * tt;
2917 tt = Lab0[2] - Lab1[2];
2918 dout[0][2] = 1.0 * tt;
2919 dout[1][2] = -1.0 * tt;
2920 rv += tt * tt;
2921 return sqrt(rv);
2922 }
2923
2924 /* Return the CIE94 Delta E color difference measure */
2925 /* including partial derivatives. */
icxdCIE94(double dout[2][3],double Lab0[3],double Lab1[3])2926 double icxdCIE94(double dout[2][3], double Lab0[3], double Lab1[3]) {
2927 double desq, _desq[2][3];
2928 double dlsq;
2929 double dcsq, _dcsq[2][2]; /* == [x][1,2] */
2930 double c12, _c12[2][2]; /* == [x][1,2] */
2931 double dhsq, _dhsq[2][2]; /* == [x][1,2] */
2932 double rv;
2933
2934 {
2935 double dl, da, db;
2936
2937 dl = Lab0[0] - Lab1[0];
2938 dlsq = dl * dl; /* dl squared */
2939 da = Lab0[1] - Lab1[1];
2940 db = Lab0[2] - Lab1[2];
2941
2942
2943 /* Compute normal Lab delta E squared */
2944 desq = dlsq + da * da + db * db;
2945 _desq[0][0] = 1.0 * dl;
2946 _desq[1][0] = -1.0 * dl;
2947 _desq[0][1] = 1.0 * da;
2948 _desq[1][1] = -1.0 * da;
2949 _desq[0][2] = 1.0 * db;
2950 _desq[1][2] = -1.0 * db;
2951 }
2952
2953 {
2954 double c1, c2, dc, tt;
2955
2956 /* Compute chromanance for the two colors */
2957 c1 = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]);
2958 c2 = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]);
2959 c12 = sqrt(c1 * c2); /* Symetric chromanance */
2960
2961 tt = 0.5 * (pow(c2, 0.5) + 1e-12)/(pow(c1, 1.5) + 1e-12);
2962 _c12[0][0] = Lab0[1] * tt;
2963 _c12[0][1] = Lab0[2] * tt;
2964 tt = 0.5 * (pow(c1, 0.5) + 1e-12)/(pow(c2, 1.5) + 1e-12);
2965 _c12[1][0] = Lab1[1] * tt;
2966 _c12[1][1] = Lab1[2] * tt;
2967
2968 /* delta chromanance squared */
2969 dc = c2 - c1;
2970 dcsq = dc * dc;
2971 if (c1 < 1e-12 || c2 < 1e-12) {
2972 c1 += 1e-12;
2973 c2 += 1e-12;
2974 }
2975 _dcsq[0][0] = -1.0 * Lab0[1] * (c2 - c1)/c1;
2976 _dcsq[0][1] = -1.0 * Lab0[2] * (c2 - c1)/c1;
2977 _dcsq[1][0] = 1.0 * Lab1[1] * (c2 - c1)/c2;
2978 _dcsq[1][1] = 1.0 * Lab1[2] * (c2 - c1)/c2;
2979 }
2980
2981 /* Compute delta hue squared */
2982 dhsq = desq - dlsq - dcsq;
2983 if (dhsq >= 0.0) {
2984 _dhsq[0][0] = _desq[0][1] - _dcsq[0][0];
2985 _dhsq[0][1] = _desq[0][2] - _dcsq[0][1];
2986 _dhsq[1][0] = _desq[1][1] - _dcsq[1][0];
2987 _dhsq[1][1] = _desq[1][2] - _dcsq[1][1];
2988 } else {
2989 dhsq = 0.0;
2990 _dhsq[0][0] = 0.0;
2991 _dhsq[0][1] = 0.0;
2992 _dhsq[1][0] = 0.0;
2993 _dhsq[1][1] = 0.0;
2994 }
2995
2996 {
2997 double sc, scsq, scf;
2998 double sh, shsq, shf;
2999
3000 /* Weighting factors for delta chromanance & delta hue */
3001 sc = 1.0 + 0.048 * c12;
3002 scsq = sc * sc;
3003
3004 sh = 1.0 + 0.014 * c12;
3005 shsq = sh * sh;
3006
3007 rv = dlsq + dcsq/scsq + dhsq/shsq;
3008
3009 scf = 0.048 * -1.0 * dcsq/(scsq * sc);
3010 shf = 0.014 * -1.0 * dhsq/(shsq * sh);
3011 dout[0][0] = _desq[0][0];
3012 dout[0][1] = _dcsq[0][0]/scsq + _c12[0][0] * scf
3013 + _dhsq[0][0]/shsq + _c12[0][0] * shf;
3014 dout[0][2] = _dcsq[0][1]/scsq + _c12[0][1] * scf
3015 + _dhsq[0][1]/shsq + _c12[0][1] * shf;
3016 dout[1][0] = _desq[1][0];
3017 dout[1][1] = _dcsq[1][0]/scsq + _c12[1][0] * scf
3018 + _dhsq[1][0]/shsq + _c12[1][0] * shf;
3019 dout[1][2] = _dcsq[1][1]/scsq + _c12[1][1] * scf
3020 + _dhsq[1][1]/shsq + _c12[1][1] * shf;
3021 return sqrt(rv);
3022 }
3023 }
3024
3025 /* ------------------------------------------------------ */
3026 /* A power-like function, based on Graphics Gems adjustment curve. */
3027 /* Avoids "toe" problem of pure power. */
3028 /* Adjusted so that "power" 2 and 0.5 agree with real power at 0.5 */
3029
icx_powlike(double vv,double pp)3030 double icx_powlike(double vv, double pp) {
3031 double tt, g;
3032
3033 if (pp >= 1.0) {
3034 g = 2.0 * (pp - 1.0);
3035 vv = vv/(g - g * vv + 1.0);
3036 } else {
3037 g = 2.0 - 2.0/pp;
3038 vv = (vv - g * vv)/(1.0 - g * vv);
3039 }
3040
3041 return vv;
3042 }
3043
3044 /* Compute the necessary aproximate power, to transform */
3045 /* the given value from src to dst. They are assumed to be */
3046 /* in the range 0.0 .. 1.0 */
icx_powlike_needed(double src,double dst)3047 double icx_powlike_needed(double src, double dst) {
3048 double pp, g;
3049
3050 if (dst <= src) {
3051 g = -((src - dst)/(dst * src - dst));
3052 pp = (0.5 * g) + 1.0;
3053 } else {
3054 g = -((src - dst)/((dst - 1.0) * src));
3055 pp = 1.0/(1.0 - 0.5 * g);
3056 }
3057
3058 return pp;
3059 }
3060
3061 /* ------------------------------------------------------ */
3062 /* Parameterized transfer/dot gain function. */
3063 /* Used for device modelling. Including partial */
3064 /* derivative for input and parameters. */
3065
3066 /* NOTE that clamping the input values seems to cause */
3067 /* conjgrad() problems. */
3068
3069 /* Transfer function */
icxTransFunc(double * v,int luord,double vv)3070 double icxTransFunc(
3071 double *v, /* Pointer to first parameter */
3072 int luord, /* Number of parameters */
3073 double vv /* Source of value */
3074 ) {
3075 double g;
3076 int ord;
3077
3078 /* Process all the shaper orders from low to high. */
3079 /* [These shapers were inspired by a Graphics Gem idea */
3080 /* (Gems IV, VI.3, "Fast Alternatives to Perlin's Bias and */
3081 /* Gain Functions, pp 401). */
3082 /* They have the nice properties that they are smooth, and */
3083 /* are monotonic. The control parameter has been */
3084 /* altered to have a range from -oo to +oo rather than 0.0 to 1.0 */
3085 /* so that the search space is less non-linear. */
3086 for (ord = 0; ord < luord; ord++) {
3087 int nsec; /* Number of sections */
3088 double sec; /* Section */
3089
3090 g = v[ord]; /* Parameter */
3091
3092 nsec = ord + 1; /* Increase sections for each order */
3093
3094 vv *= (double)nsec;
3095
3096 sec = floor(vv);
3097 if (((int)sec) & 1)
3098 g = -g; /* Alternate action in each section */
3099 vv -= sec;
3100 if (g >= 0.0) {
3101 vv = vv/(g - g * vv + 1.0);
3102 } else {
3103 vv = (vv - g * vv)/(1.0 - g * vv);
3104 }
3105 vv += sec;
3106 vv /= (double)nsec;
3107 }
3108
3109 return vv;
3110 }
3111
3112 /* Inverse transfer function */
icxInvTransFunc(double * v,int luord,double vv)3113 double icxInvTransFunc(
3114 double *v, /* Pointer to first parameter */
3115 int luord, /* Number of parameters */
3116 double vv /* Source of value */
3117 ) {
3118 double g;
3119 int ord;
3120
3121 /* Process the shaper orders in reverse from high to low. */
3122 /* [These shapers were inspired by a Graphics Gem idea */
3123 /* (Gems IV, VI.3, "Fast Alternatives to Perlin's Bias and */
3124 /* Gain Functions, pp 401). */
3125 /* They have the nice properties that they are smooth, and */
3126 /* are monotonic. The control parameter has been */
3127 /* altered to have a range from -oo to +oo rather than 0.0 to 1.0 */
3128 /* so that the search space is less non-linear. */
3129 for (ord = luord-1; ord >= 0; ord--) {
3130 int nsec; /* Number of sections */
3131 double sec; /* Section */
3132
3133 g = -v[ord]; /* Inverse parameter */
3134
3135 nsec = ord + 1; /* Increase sections for each order */
3136
3137 vv *= (double)nsec;
3138
3139 sec = floor(vv);
3140 if (((int)sec) & 1)
3141 g = -g; /* Alternate action in each section */
3142 vv -= sec;
3143 if (g >= 0.0) {
3144 vv = vv/(g - g * vv + 1.0);
3145 } else {
3146 vv = (vv - g * vv)/(1.0 - g * vv);
3147 }
3148 vv += sec;
3149 vv /= (double)nsec;
3150 }
3151
3152 return vv;
3153 }
3154
3155 /* Transfer function with offset and scale */
icxSTransFunc(double * v,int luord,double vv,double min,double max)3156 double icxSTransFunc(
3157 double *v, /* Pointer to first parameter */
3158 int luord, /* Number of parameters */
3159 double vv, /* Source of value */
3160 double min, /* Scale values */
3161 double max
3162 ) {
3163 max -= min;
3164
3165 vv = (vv - min)/max;
3166 vv = icxTransFunc(v, luord, vv);
3167 vv = (vv * max) + min;
3168 return vv;
3169 }
3170
3171 /* Inverse Transfer function with offset and scale */
icxInvSTransFunc(double * v,int luord,double vv,double min,double max)3172 double icxInvSTransFunc(
3173 double *v, /* Pointer to first parameter */
3174 int luord, /* Number of parameters */
3175 double vv, /* Source of value */
3176 double min, /* Scale values */
3177 double max
3178 ) {
3179 max -= min;
3180
3181 vv = (vv - min)/max;
3182 vv = icxInvTransFunc(v, luord, vv);
3183 vv = (vv * max) + min;
3184 return vv;
3185 }
3186
3187 /* Transfer function with partial derivative */
3188 /* with respect to the parameters. */
icxdpTransFunc(double * v,double * dv,int luord,double vv)3189 double icxdpTransFunc(
3190 double *v, /* Pointer to first parameter */
3191 double *dv, /* Return derivative wrt each parameter */
3192 int luord, /* Number of parameters */
3193 double vv /* Source of value */
3194 ) {
3195 double g;
3196 int i, ord;
3197
3198 /* Process all the shaper orders from high to low. */
3199 for (ord = 0; ord < luord; ord++) {
3200 double dsv; /* del for del in g */
3201 double ddv; /* del for del in vv */
3202 int nsec; /* Number of sections */
3203 double sec; /* Section */
3204
3205 g = v[ord]; /* Parameter */
3206
3207 nsec = ord + 1; /* Increase sections for each order */
3208
3209 vv *= (double)nsec;
3210
3211 sec = floor(vv);
3212 if (((int)sec) & 1) {
3213 g = -g; /* Alternate action in each section */
3214 }
3215 vv -= sec;
3216 if (g >= 0.0) {
3217 double tt = g - g * vv + 1.0;
3218 dsv = (vv * vv - vv)/(tt * tt);
3219 ddv = (g + 1.0)/(tt * tt);
3220 vv = vv/tt;
3221 } else {
3222 double tt = 1.0 - g * vv;
3223 dsv = (vv * vv - vv)/(tt * tt);
3224 ddv = (1.0 - g)/(tt * tt);
3225 vv = (vv - g * vv)/tt;
3226 }
3227
3228 vv += sec;
3229 vv /= (double)nsec;
3230 dsv /= (double)nsec;
3231 if (((int)sec) & 1)
3232 dsv = -dsv;
3233
3234 dv[ord] = dsv;
3235 for (i = ord - 1; i >= 0; i--)
3236 dv[i] *= ddv;
3237 }
3238
3239 return vv;
3240 }
3241
3242 /* Transfer function with offset and scale, and */
3243 /* partial derivative with respect to the parameters. */
icxdpSTransFunc(double * v,double * dv,int luord,double vv,double min,double max)3244 double icxdpSTransFunc(
3245 double *v, /* Pointer to first parameter */
3246 double *dv, /* Return derivative wrt each parameter */
3247 int luord, /* Number of parameters */
3248 double vv, /* Source of value */
3249 double min, /* Scale values */
3250 double max
3251 ) {
3252 int i;
3253 max -= min;
3254
3255 vv = (vv - min)/max;
3256 vv = icxdpTransFunc(v, dv, luord, vv);
3257 vv = (vv * max) + min;
3258 for (i = 0; i < luord; i++)
3259 dv[i] *= max;
3260 return vv;
3261 }
3262
3263
3264 /* Transfer function with partial derivative */
3265 /* with respect to the input value. */
icxdiTransFunc(double * v,double * pdin,int luord,double vv)3266 double icxdiTransFunc(
3267 double *v, /* Pointer to first parameter */
3268 double *pdin, /* Return derivative wrt source value */
3269 int luord, /* Number of parameters */
3270 double vv /* Source of value */
3271 ) {
3272 double g, din;
3273 int ord;
3274
3275 #ifdef NEVER
3276 if (vv < 0.0 || vv > 1.0) {
3277 if (vv < 0.0)
3278 vv = 0.0;
3279 else
3280 vv = 1.0;
3281
3282 *pdin = 0.0;
3283 return vv;
3284 }
3285 #endif
3286 din = 1.0;
3287
3288 /* Process all the shaper orders from high to low. */
3289 for (ord = 0; ord < luord; ord++) {
3290 double ddv; /* del for del in vv */
3291 int nsec; /* Number of sections */
3292 double sec; /* Section */
3293
3294 g = v[ord]; /* Parameter */
3295
3296 nsec = ord + 1; /* Increase sections for each order */
3297
3298 vv *= (double)nsec;
3299
3300 sec = floor(vv);
3301 if (((int)sec) & 1) {
3302 g = -g; /* Alternate action in each section */
3303 }
3304 vv -= sec;
3305 if (g >= 0.0) {
3306 double tt = g - g * vv + 1.0;
3307 ddv = (g + 1.0)/(tt * tt);
3308 vv = vv/tt;
3309 } else {
3310 double tt = 1.0 - g * vv;
3311 ddv = (1.0 - g)/(tt * tt);
3312 vv = (vv - g * vv)/tt;
3313 }
3314
3315 vv += sec;
3316 vv /= (double)nsec;
3317 din *= ddv;
3318 }
3319
3320 *pdin = din;
3321 return vv;
3322 }
3323
3324 /* Transfer function with offset and scale, and */
3325 /* partial derivative with respect to the input value. */
icxdiSTransFunc(double * v,double * pdv,int luord,double vv,double min,double max)3326 double icxdiSTransFunc(
3327 double *v, /* Pointer to first parameter */
3328 double *pdv, /* Return derivative wrt source value */
3329 int luord, /* Number of parameters */
3330 double vv, /* Source of value */
3331 double min, /* Scale values */
3332 double max
3333 ) {
3334 max -= min;
3335
3336 vv = (vv - min)/max;
3337 vv = icxdiTransFunc(v, pdv, luord, vv);
3338 vv = (vv * max) + min;
3339 return vv;
3340 }
3341
3342
3343 /* Transfer function with partial derivative */
3344 /* with respect to the parameters and the input value. */
icxdpdiTransFunc(double * v,double * dv,double * pdin,int luord,double vv)3345 double icxdpdiTransFunc(
3346 double *v, /* Pointer to first parameter */
3347 double *dv, /* Return derivative wrt each parameter */
3348 double *pdin, /* Return derivative wrt source value */
3349 int luord, /* Number of parameters */
3350 double vv /* Source of value */
3351 ) {
3352 double g, din;
3353 int i, ord;
3354
3355 #ifdef NEVER
3356 if (vv < 0.0 || vv > 1.0) {
3357 if (vv < 0.0)
3358 vv = 0.0;
3359 else
3360 vv = 1.0;
3361
3362 for (ord = 0; ord < luord; ord++)
3363 dv[ord] = 0.0;
3364 *pdin = 0.0;
3365 return vv;
3366 }
3367 #endif
3368 din = 1.0;
3369
3370 /* Process all the shaper orders from high to low. */
3371 for (ord = 0; ord < luord; ord++) {
3372 double dsv; /* del for del in g */
3373 double ddv; /* del for del in vv */
3374 int nsec; /* Number of sections */
3375 double sec; /* Section */
3376
3377 g = v[ord]; /* Parameter */
3378
3379 nsec = ord + 1; /* Increase sections for each order */
3380
3381 vv *= (double)nsec;
3382
3383 sec = floor(vv);
3384 if (((int)sec) & 1) {
3385 g = -g; /* Alternate action in each section */
3386 }
3387 vv -= sec;
3388 if (g >= 0.0) {
3389 double tt = g - g * vv + 1.0;
3390 dsv = (vv * vv - vv)/(tt * tt);
3391 ddv = (g + 1.0)/(tt * tt);
3392 vv = vv/tt;
3393 } else {
3394 double tt = 1.0 - g * vv;
3395 dsv = (vv * vv - vv)/(tt * tt);
3396 ddv = (1.0 - g)/(tt * tt);
3397 vv = (vv - g * vv)/tt;
3398 }
3399
3400 vv += sec;
3401 vv /= (double)nsec;
3402 dsv /= (double)nsec;
3403 if (((int)sec) & 1)
3404 dsv = -dsv;
3405
3406 dv[ord] = dsv;
3407 for (i = ord - 1; i >= 0; i--)
3408 dv[i] *= ddv;
3409 din *= ddv;
3410 }
3411
3412 *pdin = din;
3413 return vv;
3414 }
3415
3416 /* Transfer function with offset and scale, and */
3417 /* partial derivative with respect to the */
3418 /* parameters and the input value. */
icxdpdiSTransFunc(double * v,double * dv,double * pdin,int luord,double vv,double min,double max)3419 double icxdpdiSTransFunc(
3420 double *v, /* Pointer to first parameter */
3421 double *dv, /* Return derivative wrt each parameter */
3422 double *pdin, /* Return derivative wrt source value */
3423 int luord, /* Number of parameters */
3424 double vv, /* Source of value */
3425 double min, /* Scale values */
3426 double max
3427 ) {
3428 int i;
3429 max -= min;
3430
3431 vv = (vv - min)/max;
3432 vv = icxdpdiTransFunc(v, dv, pdin, luord, vv);
3433 vv = (vv * max) + min;
3434 for (i = 0; i < luord; i++)
3435 dv[i] *= max;
3436 return vv;
3437 }
3438
3439 /* ------------------------------------------------------ */
3440 /* Multi-plane interpolation, used for device modelling. */
3441 /* Including partial derivative for input and parameters. */
3442 /* A simple flat plane is used for each output. */
3443
3444 /* Multi-plane interpolation - uses base + di slope values. */
3445 /* Parameters are assumed to be fdi groups of di + 1 parameters. */
icxPlaneInterp(double * v,int fdi,int di,double * out,double * in)3446 void icxPlaneInterp(
3447 double *v, /* Pointer to first parameter [fdi * (di + 1)] */
3448 int fdi, /* Number of output channels */
3449 int di, /* Number of input channels */
3450 double *out, /* Resulting fdi values */
3451 double *in /* Input di values */
3452 ) {
3453 int e, f;
3454
3455 for (f = 0; f < fdi; f++) {
3456 for (out[f] = 0.0, e = 0; e < di; e++, v++) {
3457 out[f] += in[e] * *v;
3458 }
3459 out[f] += *v;
3460 }
3461 }
3462
3463
3464 /* Multii-plane interpolation with partial derivative */
3465 /* with respect to the input and parameters. */
icxdpdiPlaneInterp(double * v,double * dv,double * din,int fdi,int di,double * out,double * in)3466 void icxdpdiPlaneInterp(
3467 double *v, /* Pointer to first parameter value [fdi * (di + 1)] */
3468 double *dv, /* Return [1 + di] deriv. wrt each parameter v */
3469 double *din, /* Return [fdi * di] deriv. wrt each input value */
3470 int fdi, /* Number of output channels */
3471 int di, /* Number of input channels */
3472 double *out, /* Resulting fdi values */
3473 double *in /* Input di values */
3474 ) {
3475 int e, ee, f, g;
3476 int dip2 = (di + 1); /* Output dim increment through parameters */
3477
3478 /* Compute the output values */
3479 for (f = 0; f < fdi; f++) {
3480 for (out[f] = 0.0, e = 0; e < di; e++)
3481 out[f] += in[e] * v[f * dip2 + e];
3482 out[f] += v[f * dip2 + e];
3483 }
3484
3485 /* Since interpolation is verys simple, derivative are also simple */
3486
3487 /* Copy del for parameter to return array */
3488 for (e = 0; e < di; e++)
3489 dv[e] = in[e];
3490 dv[e] = 1.0;
3491
3492 /* Compute del of out[] from in[] */
3493 for (f = 0; f < fdi; f++) {
3494 for (e = 0; e < di; e++) {
3495 din[f * di + e] = v[f * dip2 + e];
3496 }
3497 }
3498 }
3499
3500 /* ------------------------------------------------------ */
3501 /* Matrix cube interpolation, used for device modelling. */
3502 /* Including partial derivative for input and parameters. */
3503
3504 /* Matrix cube interpolation - interpolate between 2^di output corner values. */
3505 /* Parameters are assumed to be fdi groups of 2^di parameters. */
icxCubeInterp(double * v,int fdi,int di,double * out,double * in)3506 void icxCubeInterp(
3507 double *v, /* Pointer to first parameter */
3508 int fdi, /* Number of output channels */
3509 int di, /* Number of input channels */
3510 double *out, /* Resulting fdi values */
3511 double *in /* Input di values */
3512 ) {
3513 int e, f, g;
3514 double gw[1 << MXDI]; /* weight for each matrix grid cube corner */
3515
3516 /* Compute corner weights needed for interpolation */
3517 gw[0] = 1.0;
3518 for (e = 0, g = 1; e < di; e++, g *= 2) {
3519 int i;
3520 for (i = 0; i < g; i++) {
3521 gw[g+i] = gw[i] * in[e];
3522 gw[i] *= (1.0 - in[e]);
3523 }
3524 }
3525
3526 /* Now compute the output values */
3527 for (f = 0; f < fdi; f++) {
3528 out[f] = 0.0; /* For each output value */
3529 for (e = 0; e < (1 << di); e++) { /* For all corners of cube */
3530 out[f] += gw[e] * *v;
3531 v++;
3532 }
3533 }
3534 }
3535
3536
3537 /* Matrix cube interpolation. with partial derivative */
3538 /* with respect to the input and parameters. */
icxdpdiCubeInterp(double * v,double * dv,double * din,int fdi,int di,double * out,double * in)3539 void icxdpdiCubeInterp(
3540 double *v, /* Pointer to first parameter value [fdi * 2^di] */
3541 double *dv, /* Return [2^di] deriv. wrt each parameter v */
3542 double *din, /* Return [fdi * di] deriv. wrt each input value */
3543 int fdi, /* Number of output channels */
3544 int di, /* Number of input channels */
3545 double *out, /* Resulting fdi values */
3546 double *in /* Input di values */
3547 ) {
3548 int e, ee, f, g;
3549 int dip2 = (1 << di);
3550 double gw[1 << MXDI]; /* weight for each matrix grid cube corner */
3551
3552 /* Compute corner weights needed for interpolation */
3553 gw[0] = 1.0;
3554 for (e = 0, g = 1; e < di; e++, g *= 2) {
3555 int i;
3556 for (i = 0; i < g; i++) {
3557 gw[g+i] = gw[i] * in[e];
3558 gw[i] *= (1.0 - in[e]);
3559 }
3560 }
3561
3562 /* Now compute the output values */
3563 for (f = 0; f < fdi; f++) {
3564 out[f] = 0.0; /* For each output value */
3565 for (ee = 0; ee < dip2; ee++) { /* For all corners of cube */
3566 out[f] += gw[ee] * v[f * dip2 + ee];
3567 }
3568 }
3569
3570 /* Copy del for parameter to return array */
3571 for (ee = 0; ee < dip2; ee++) { /* For all other corners of cube */
3572 dv[ee] = gw[ee]; /* del from parameter */
3573 }
3574
3575 /* Compute del from in[] value we want */
3576 for (e = 0; e < di; e++) { /* For input we want del wrt */
3577
3578 for (f = 0; f < fdi; f++)
3579 din[f * di + e] = 0.0; /* Zero del ready of accumulation */
3580
3581 for (ee = 0; ee < dip2; ee++) { /* For all corners of cube weights, */
3582 int e2; /* accumulate del from in[] we want. */
3583 double vv = 1.0;
3584
3585 /* Compute in[] weighted cube corners for all except del of in[] we want */
3586 for (e2 = 0; e2 < di; e2++) { /* const from non del inputs */
3587 if (e2 == e)
3588 continue;
3589 if (ee & (1 << e2))
3590 vv *= in[e2];
3591 else
3592 vv *= (1.0 - in[e2]);
3593 }
3594
3595 /* Accumulate contribution of in[] we want for corner to out[] we want */
3596 if (ee & (1 << e)) {
3597 for (f = 0; f < fdi; f++)
3598 din[f * di + e] += v[f * dip2 + ee] * vv;
3599 } else {
3600 for (f = 0; f < fdi; f++)
3601 din[f * di + e] -= v[f * dip2 + ee] * vv;
3602 }
3603 }
3604 }
3605 }
3606
3607 /* ------------------------------------------------------ */
3608 /* Matrix cube simplex interpolation, used for device modelling. */
3609
3610 /* Matrix cube simplex interpolation - interpolate between 2^di output corner values. */
3611 /* Parameters are assumed to be fdi groups of 2^di parameters. */
icxCubeSxInterp(double * v,int fdi,int di,double * out,double * in)3612 void icxCubeSxInterp(
3613 double *v, /* Pointer to first parameter */
3614 int fdi, /* Number of output channels */
3615 int di, /* Number of input channels */
3616 double *out, /* Resulting fdi values */
3617 double *in /* Input di values */
3618 ) {
3619 int si[MAX_CHAN]; /* in[] Sort index, [0] = smalest */
3620
3621 //{
3622 // double tout[MXDO];
3623 //
3624 // icxCubeInterp(v, fdi, di, tout, in);
3625 //printf("\n~1 Cube interp result = %f\n",tout[0]);
3626 //}
3627
3628 //printf("~1 icxCubeSxInterp: %f %f %f\n", in[0], in[1], in[2]);
3629 /* Do insertion sort on coordinates, smallest to largest. */
3630 {
3631 int ff, vf;
3632 unsigned int e;
3633 double v;
3634 for (e = 0; e < di; e++)
3635 si[e] = e; /* Initial unsorted indexes */
3636
3637 for (e = 1; e < di; e++) {
3638 ff = e;
3639 v = in[si[ff]];
3640 vf = ff;
3641 while (ff > 0 && in[si[ff-1]] > v) {
3642 si[ff] = si[ff-1];
3643 ff--;
3644 }
3645 si[ff] = vf;
3646 }
3647 }
3648 //printf("~1 sort order %d %d %d\n", si[0], si[1], si[2]);
3649 //printf(" from %f %f %f\n", in[si[0]], in[si[1]], in[si[2]]);
3650
3651 /* Now compute the weightings, simplex vertices and output values */
3652 {
3653 unsigned int e, f;
3654 double w; /* Current vertex weight */
3655
3656 w = 1.0 - in[si[di-1]]; /* Vertex at base of cell */
3657 for (f = 0; f < fdi; f++) {
3658 out[f] = w * v[f * (1 << di)];
3659 //printf("~1 out[%d] = %f = %f * %f\n",f,out[f],w,v[f * (1 << di)]);
3660 }
3661
3662 for (e = di-1; e > 0; e--) { /* Middle verticies */
3663 w = in[si[e]] - in[si[e-1]];
3664 v += (1 << si[e]); /* Move to top of cell in next largest dimension */
3665 for (f = 0; f < fdi; f++) {
3666 out[f] += w * v[f * (1 << di)];
3667 //printf("~1 out[%d] = %f += %f * %f\n",f,out[f],w,v[f * (1 << di)]);
3668 }
3669 }
3670
3671 w = in[si[0]];
3672 v += (1 << si[0]); /* Far corner from base of cell */
3673 for (f = 0; f < fdi; f++) {
3674 out[f] += w * v[f * (1 << di)];
3675 //printf("~1 out[%d] = %f += %f * %f\n",f,out[f],w,v[f * (1 << di)]);
3676 }
3677 }
3678 }
3679
3680 /* ------------------------------------------------------ */
3681 /* Matrix multiplication, used for device modelling. */
3682 /* Including partial derivative for input and parameters. */
3683
3684
3685 /* 3x3 matrix in 1D array multiplication */
icxMulBy3x3Parm(double out[3],double mat[9],double in[3])3686 void icxMulBy3x3Parm(
3687 double out[3], /* Return input multiplied by matrix */
3688 double mat[9], /* Matrix organised in [slow][fast] order */
3689 double in[3] /* Input values */
3690 ) {
3691 double *v = mat, ov[3];
3692 int e, f;
3693
3694 /* Compute the output values */
3695 for (f = 0; f < 3; f++) {
3696 ov[f] = 0.0; /* For each output value */
3697 for (e = 0; e < 3; e++) {
3698 ov[f] += *v++ * in[e];
3699 }
3700 }
3701 out[0] = ov[0];
3702 out[1] = ov[1];
3703 out[2] = ov[2];
3704 }
3705
3706
3707 /* 3x3 matrix in 1D array multiplication, with partial derivatives */
3708 /* with respect to just the input. */
icxdpdiiMulBy3x3Parm(double out[3],double din[3][3],double mat[9],double in[3])3709 void icxdpdiiMulBy3x3Parm(
3710 double out[3], /* Return input multiplied by matrix */
3711 double din[3][3], /* Return deriv for each [output] with respect to [input] */
3712 double mat[9], /* Matrix organised in [slow][fast] order */
3713 double in[3] /* Input values */
3714 ) {
3715 double *v, ov[3];
3716 int e, f;
3717
3718 /* Compute the output values */
3719 v = mat;
3720 for (f = 0; f < 3; f++) {
3721 ov[f] = 0.0; /* For each output value */
3722 for (e = 0; e < 3; e++) {
3723 ov[f] += *v++ * in[e];
3724 }
3725 }
3726
3727 /* Compute deriv. with respect to the input values */
3728 /* This is pretty simple for a matrix ... */
3729 v = mat;
3730 for (f = 0; f < 3; f++)
3731 for (e = 0; e < 3; e++)
3732 din[f][e] = *v++;
3733
3734 out[0] = ov[0];
3735 out[1] = ov[1];
3736 out[2] = ov[2];
3737 }
3738
3739 /* 3x3 matrix in 1D array multiplication, with partial derivatives */
3740 /* with respect to the input and parameters. */
icxdpdiMulBy3x3Parm(double out[3],double dv[3][9],double din[3][3],double mat[9],double in[3])3741 void icxdpdiMulBy3x3Parm(
3742 double out[3], /* Return input multiplied by matrix */
3743 double dv[3][9], /* Return deriv for each [output] with respect to [param] */
3744 double din[3][3], /* Return deriv for each [output] with respect to [input] */
3745 double mat[9], /* Matrix organised in [slow][fast] order */
3746 double in[3] /* Input values */
3747 ) {
3748 double *v, ov[3];
3749 int e, f;
3750
3751 /* Compute the output values */
3752 v = mat;
3753 for (f = 0; f < 3; f++) {
3754 ov[f] = 0.0; /* For each output value */
3755 for (e = 0; e < 3; e++) {
3756 ov[f] += *v++ * in[e];
3757 }
3758 }
3759
3760 /* Compute deriv. with respect to the matrix parameter % 3 */
3761 /* This is pretty simple for a matrix ... */
3762 for (f = 0; f < 3; f++) {
3763 for (e = 0; e < 9; e++) {
3764 if (e/3 == f)
3765 dv[f][e] = in[e % 3];
3766 else
3767 dv[f][e] = 0.0;
3768 }
3769 }
3770
3771 /* Compute deriv. with respect to the input values */
3772 /* This is pretty simple for a matrix ... */
3773 v = mat;
3774 for (f = 0; f < 3; f++)
3775 for (e = 0; e < 3; e++)
3776 din[f][e] = *v++;
3777
3778 out[0] = ov[0];
3779 out[1] = ov[1];
3780 out[2] = ov[2];
3781 }
3782
3783 /* - - - - - - - - - - */
3784 #undef stricmp
3785
3786
3787
3788
3789
3790
3791
3792
3793
3794
3795
3796
3797
3798
3799
3800
3801
3802
3803
3804
3805
3806