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 * This is the third major version of xlut.c (originally called xlut2.c)
15 * Based on the old xlut.c code (preserved as xlut1.c)
16 * This version uses xfit.c to do the curve and rspl fitting.
17 */
18
19 /*
20 * This module provides the expanded icclib functionality
21 * for lut based profiles.
22 * This file is #included in xicc.c, to keep its functions private.
23 *
24 * This version creates both input and output 1D luts by
25 * optimising the accuracy of the profile for a linear clut.
26 *
27 */
28
29 /*
30 TTBD:
31
32 See gamut/gammap.c for notes on the desirability of
33 a non-minimum delta E colorimetric intent as the default.
34
35 Should fix gamut's so that they have enough information
36 to spefify a gamut mapping without the need for
37 a source colorspace profile, and fix the code
38 to work with just a gamut. This would take us further
39 towards supporting the PRMG reference gamut interoperability
40 as an option.
41
42 Should the input profile white point determination
43 be made a bit smarter about determining the chromaticity ?
44 ie. not take on the tint of the whitest patch, but an
45 average of neutral patches ?
46
47 Should xlutfix.c be revived (also adding ICM_CLUT_SET_APXLS support),
48 to improve "bumpy black" problem ?
49
50 Would be nice to be able to specify a specific patch
51 as the white one rather than using heuristic to identify it,
52 since some pathalogical cases don't work.
53 */
54
55 /*
56
57 !!!!! This has all been fixed ? !!!!!!
58
59 NOTE :- an alternative to the way display profile absolute is handled here
60 would be to always chromatically adapt the illuminant to D50, and encode
61 that in the Chromatic adapation tag. To make absolute colorimetric
62 do anything useful though, the chromatic adapation tag would
63 have to be used for absolute intent.
64 This may be the way of improving compatibility with other systems,
65 and is needed for V4, but would break compatibility with existing
66 Argyll profiles, unless special measures are taken:
67
68 ie.
69
70 1) if (display profile & using chromatic adaptation tag)
71 Create Bradford chromatic adapation matrix and store it in tag
72 Adapt all the readings using Bradford
73 Create white point and store it in tag (white point will be D50)
74 Adapt all the readings to the white point using wrong Von-Kries (== NOOP)
75 Store relative colorimetric cLUT
76 Set version >= 2.4
77
78 else
79 2) if (display scheme A or using Argyll historical printer scheme)
80 Create white point and store it in tag
81 Adapt all the readings to the white point using Bradford
82 Store relative colorimetric tag
83 Set version < 2.4 for V2 profile
84 Add private Absolute Transform Matrix (labels V4 profile)
85
86 3) else (display scheme B or strict ICC printer compatibility)
87 Create white point and store it in tag
88 Adapt all the readings to the white point using Wrong Von-Kries
89 Store relative colorimetric tag
90 Set version >= 2.4
91
92
93 Argyll Processing for each type
94
95 1) if display and chromatic adapation matrix
96 Un-adapt matrix or cLUT using wrong Von-Kries from white point
97 Un-adapt matrix or cLUT using chromatic matrix
98 Un-adapt apparant white point & black point using chromatic transform
99
100 if (not absolute intent)
101 Create Bradford transfor from white to PCS D50
102 Adapt all matrix or cLUT
103 else
104 2) if (display scheme A or using Argyll < V2.4. profile
105 or find Absolute Transform Matrix)
106 if (absolute intent)
107 Un-adapt matrix or cLUT using Bradford from white point
108
109 3) else (display scheme B or !Argyll profile or ( >= V2.4 profile
110 and !Absolute Transform Matrix))
111 Un-adapt matrix or cLUT using wrong Von-Kries from white point
112
113 if (not absolute intent)
114 Create Bradford transfor from white to PCS D50
115 Adapt all matrix or cLUT to white
116
117
118 The problem with this is that it wouldn't do the right thing on old Argyll
119 type profiles that weren't labeled or recognized.
120
121 Is there a way of recognizing Bradford Absolute transform Matricies if
122 the color chromaticities are given ?
123
124 */
125
126 /*
127 A similar condrum is that it seems that an unwritten convention for
128 V2 profiles is to scale the black point of the perceptual and
129 saturation tables to 0 (Part of the V4 spec is to scale to Y = 3.1373).
130
131 To get better gamut mapping we should therefore unscale the perceptual
132 and saturation A2B table to have the same black point as the colorimetric
133 table before computing the gamut mapping, and then apply the opposite
134 transform to our perceptual B2A and A2B tables.
135
136 */
137
138 /*
139 There is interesting behaviour for Jab 0,0,0 in, in that
140 it gets mapped to (something like) Lab -1899.019855 -213.574625 -231.914285
141 which after per-component clipping of the inv out tables looks like
142 0, -128, -128, which may be mapped to (say) RGB 0.085455 1.000000 0.936951,
143 which is not black.
144
145 */
146
147 #include "xfit.h"
148
149 #undef USE_CIE94_DE /* [Undef] Use CIE94 delta E measure when creating in/out curves */
150 /* Don't use CIE94 because it makes peak error worse ? */
151
152 #undef DEBUG /* [Undef] Verbose debug information */
153 #undef CNDTRACE /* [Undef] Enable xicc->trace conditional verbosity */
154 #undef DEBUG_OLUT /* [Undef] Print steps in overall fwd & rev lookups */
155 #undef DEBUG_RLUT /* [Undef] Print values being reverse lookup up */
156 #undef DEBUG_SPEC /* [Undef] Debug some specific cases */
157 #undef INK_LIMIT_TEST /* [Undef] Turn off input tables for ink limit testing */
158 #undef CHECK_ILIMIT /* [Undef] Do sanity checks on meeting ink limit */
159 #undef WARN_CLUT_CLIPPING /* [Undef] Print warning if setting clut clips */
160 #undef DISABLE_KCURVE_FILTER /* [Undef] don't filter the Kcurve */
161 #undef REPORT_LOCUS_SEGMENTS /* [Undef[ Examine how many segments there are in aux inversion */
162
163 #undef FASTREVSETUP_NON_CAM /* [Undef] Use fast setup on innerm non-CAM lookup, if we're */
164 /* going to use CAM clip for nn lookup */
165
166 #define XYZ_EXTRA_SMOOTH 20.0 /* Extra smoothing factor for XYZ profiles */
167 /* !!! Note this is mainly due to smoothing being */
168 /* scaled by data range in rspl code !!! */
169 #define SHP_SMOOTH 1.0 /* Input shaper curve smoothing */
170 #define OUT_SMOOTH1 1.0 /* Output shaper curve smoothing for L*, X,Y,Z */
171 #define OUT_SMOOTH2 1.0 /* Output shaper curve smoothing for a*, b* */
172
173 #define CAMCLIPTRANS 1.0 /* [1.0] Cam clipping transition region Delta E */
174 /* Should this be smaller ? */
175 #undef USECAMCLIPSPLINE /* [Und] use spline blend between PCS and Jab */
176
177 #define USELCHWEIGHT /* [def] Use LCh weighting for nn clip if possible */
178 #define JCCWEIGHT 2.0 /* [2.0] Amount to emphasize J delta E in in computing clip */
179 #define CCCWEIGHT 1.0 /* [1.0] Amount to emphasize C delta E in in computing clip */
180 #define HCCWEIGHT 2.2 /* [2.2] Amount to emphasize H delta E in in computing clip */
181
182 /*
183 * TTBD:
184 *
185 * Reverse lookup of Lab
186 * Make NEARCLIP the default ??
187 *
188 * XYZ vector clipping isn't implemented.
189 *
190 * Some of the error handling is crude. Shouldn't use
191 * error(), should return status.
192 */
193
194 static double icxLimitD(icxLuLut *p, double *in); /* For input' */
195 #define icxLimitD_void ((double (*)(void *, double *))icxLimitD) /* Cast with void 1st arg */
196 static double icxLimit(icxLuLut *p, double *in); /* For input */
197 static int icxLuLut_init_clut_camclip(icxLuLut *p);
198
199 /* Debug overall lookup */
200 #ifdef DEBUG_OLUT
201 #undef DBOL
202 #ifdef CNDTRACE
203 #define DBOL(xxx) if (p->trace) printf xxx ;
204 #else
205 #define DBOL(xxx) printf xxx ;
206 #endif
207 #else
208 #undef DBOL
209 #define DBOL(xxx)
210 #endif
211
212 /* Debug reverse lookup */
213 #ifdef DEBUG_RLUT
214 #undef DBR
215 #ifdef CNDTRACE
216 #define DBR(xxx) if (p->trace) printf xxx ;
217 #else
218 #define DBR(xxx) printf xxx ;
219 #endif
220 #else
221 #undef DBR
222 #define DBR(xxx)
223 #endif
224
225 /* Debug some specific cases (fwd_relpcs_outpcs, bwd_outpcs_relpcs) */
226 #ifdef DEBUG_SPEC
227 # undef DBS
228 # ifdef CNDTRACE
229 # define DBS(xxx) if (p->trace) printf xxx ;
230 # else
231 # define DBS(xxx) printf xxx ;
232 # endif
233 #else
234 # undef DBS
235 # define DBS(xxx)
236 #endif
237
238 /* ========================================================== */
239 /* xicc lookup code. */
240 /* ========================================================== */
241
242 /* Forward and Backward Multi-Dimensional Interpolation type conversion */
243 /* Return 0 on success, 1 if clipping occured, 2 on other error */
244
245 /* Components of overall lookup, in order */
246
icxLuLut_in_abs(icxLuLut * p,double * out,double * in)247 int icxLuLut_in_abs(icxLuLut *p, double *out, double *in) {
248 int rv = 0;
249
250 if (p->ins == icxSigJabData) {
251 DBOL(("xlut in_abs: CAM in = %s\n", icmPdv(p->inputChan, in)));
252 p->cam->cam_to_XYZ(p->cam, out, in);
253 DBOL(("xlut in_abs: XYZ = %s\n", icmPdv(p->inputChan, out)));
254 /* Hack to prevent CAM02 weirdness being amplified by inv_abs() */
255 /* or any later per channel clipping. */
256 /* Limit -Y to non-stupid values by scaling */
257 if (out[1] < -0.1) {
258 out[0] *= -0.1/out[1];
259 out[2] *= -0.1/out[1];
260 out[1] = -0.1;
261 DBOL(("xlut in_abs: after clipping -Y %s\n",icmPdv(p->outputChan, out)));
262 }
263 rv |= ((icmLuLut *)p->plu)->in_abs((icmLuLut *)p->plu, out, out);
264 DBOL(("xlut in_abs: XYZ out = %s\n", icmPdv(p->inputChan, out)));
265 } else {
266 DBOL(("xlut in_abs: PCS in = %s\n", icmPdv(p->inputChan, in)));
267 rv |= ((icmLuLut *)p->plu)->in_abs((icmLuLut *)p->plu, out, in);
268 DBOL(("xlut in_abs: PCS out = %s\n", icmPdv(p->inputChan, out)));
269 }
270
271 return rv;
272 }
273
274 /* Possible matrix lookup */
275 /* input->input (not distinguishing matrix altered input values) */
icxLuLut_matrix(icxLuLut * p,double * out,double * in)276 int icxLuLut_matrix(icxLuLut *p, double *out, double *in) {
277 int rv = 0;
278 rv |= ((icmLuLut *)p->plu)->matrix((icmLuLut *)p->plu, out, in);
279 return rv;
280 }
281
282 /* Do input -> input' lookup */
icxLuLut_input(icxLuLut * p,double * out,double * in)283 int icxLuLut_input(icxLuLut *p, double *out, double *in) {
284 #ifdef NEVER
285 return ((icmLuLut *)p->plu)->input((icmLuLut *)p->plu, out, in);
286 #else
287 int rv = 0;
288 co tc;
289 int i;
290 for (i = 0; i < p->inputChan; i++) {
291 tc.p[0] = in[i];
292 rv |= p->inputTable[i]->interp(p->inputTable[i], &tc);
293 out[i] = tc.v[0];
294 }
295 return rv;
296 #endif
297 }
298
299 /* Do input'->output' lookup, with aux' return */
300 /* (The aux' is just extracted from the in' values) */
icxLuLut_clut_aux(icxLuLut * p,double * out,double * oink,double * auxv,double * in)301 int icxLuLut_clut_aux(icxLuLut *p,
302 double *out, /* output' value */
303 double *oink, /* If not NULL, return amount input is over the ink limit, 0 if not */
304 double *auxv, /* If not NULL, return aux value used (packed) */
305 double *in /* input' value */
306 ) {
307 int rv = 0;
308 co tc;
309 int i;
310
311 for (i = 0; i < p->inputChan; i++)
312 tc.p[i] = in[i];
313 rv |= p->clutTable->interp(p->clutTable, &tc);
314 for (i = 0; i < p->outputChan; i++)
315 out[i] = tc.v[i];
316
317 if (auxv != NULL) {
318 int ee = 0;
319 for (i = 0; i < p->clutTable->di; i++) {
320 double v = in[i];
321 if (p->auxm[i] != 0) {
322 auxv[ee] = v;
323 ee++;
324 }
325 }
326 }
327
328 if (oink != NULL) {
329 double lim = 0.0;
330
331 if (p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0) {
332 lim = icxLimitD(p, in);
333 if (lim < 0.0)
334 lim = 0.0;
335 }
336 *oink = lim;
337 }
338
339 return rv;
340 }
341
342 /* Do input'->output' lookup */
icxLuLut_clut(icxLuLut * p,double * out,double * in)343 int icxLuLut_clut(icxLuLut *p, double *out, double *in) {
344 #ifdef NEVER
345 return ((icmLuLut *)p->plu)->clut((icmLuLut *)p->plu, out, in);
346 #else
347 return icxLuLut_clut_aux(p, out, NULL, NULL, in);
348 #endif
349 }
350
351 /* Do output'->output lookup */
icxLuLut_output(icxLuLut * p,double * out,double * in)352 int icxLuLut_output(icxLuLut *p, double *out, double *in) {
353 int rv = 0;
354
355 if (p->mergeclut == 0) {
356 #ifdef NEVER
357 rv = ((icmLuLut *)p->plu)->output((icmLuLut *)p->plu, out, in);
358 #else
359 co tc;
360 int i;
361 for (i = 0; i < p->outputChan; i++) {
362 tc.p[0] = in[i];
363 rv |= p->outputTable[i]->interp(p->outputTable[i], &tc);
364 out[i] = tc.v[0];
365 }
366 #endif
367 } else {
368 int i;
369 for (i = 0; i < p->outputChan; i++)
370 out[i] = in[i];
371 }
372 return rv;
373 }
374
375 /* Relative to absolute conversion + PCS to PCS override (Effective PCS) conversion */
icxLuLut_out_abs(icxLuLut * p,double * out,double * in)376 int icxLuLut_out_abs(icxLuLut *p, double *out, double *in) {
377 int rv = 0;
378 if (p->mergeclut == 0) {
379 DBOL(("xlut out_abs: PCS in = %s\n", icmPdv(p->outputChan, in)));
380
381 rv |= ((icmLuLut *)p->plu)->out_abs((icmLuLut *)p->plu, out, in);
382
383 DBOL(("xlut out_abs: ABS PCS out = %s\n", icmPdv(p->outputChan, out)));
384
385 if (p->outs == icxSigJabData) {
386 p->cam->XYZ_to_cam(p->cam, out, out);
387
388 DBOL(("xlut out_abs: CAM out = %s\n", icmPdv(p->outputChan, out)));
389 }
390 } else {
391 int i;
392 for (i = 0; i < p->outputChan; i++)
393 out[i] = in[i];
394 }
395
396 return rv;
397 }
398
399 /* Overall lookup */
400 static int
icxLuLut_lookup(icxLuBase * pp,double * out,double * in)401 icxLuLut_lookup (
402 icxLuBase *pp, /* This */
403 double *out, /* Vector of output values */
404 double *in /* Vector of input values */
405 ) {
406 icxLuLut *p = (icxLuLut *)pp;
407 int rv = 0;
408 double temp[MAX_CHAN];
409
410 DBOL(("xicclu: in = %s\n", icmPdv(p->inputChan, in)));
411 rv |= p->in_abs (p, temp, in);
412 DBOL(("xicclu: after abs = %s\n", icmPdv(p->inputChan, temp)));
413 rv |= p->matrix (p, temp, temp);
414 DBOL(("xicclu: after matrix = %s\n", icmPdv(p->inputChan, temp)));
415 rv |= p->input (p, temp, temp);
416 DBOL(("xicclu: after inout = %s\n", icmPdv(p->inputChan, temp)));
417 rv |= p->clut (p, out, temp);
418 DBOL(("xicclu: after clut = %s\n", icmPdv(p->outputChan, out)));
419 if (p->mergeclut == 0) {
420 rv |= p->output (p, out, out);
421 DBOL(("xicclu: after output = %s\n", icmPdv(p->outputChan, out)));
422 rv |= p->out_abs (p, out, out);
423 DBOL(("xicclu: after outabs = %s\n", icmPdv(p->outputChan, out)));
424 }
425 return rv;
426 }
427
428 /* - - - - - - - - - - - - - - - - - - - - - - - - - - */
429 /* Given a relative XYZ or Lab PCS value, convert in the fwd direction into */
430 /* the nominated output PCS (ie. Absolute, Jab etc.) */
431 /* (This is used in generating gamut compression in B2A tables) */
icxLuLut_fwd_relpcs_outpcs(icxLuBase * pp,icColorSpaceSignature is,double * out,double * in)432 void icxLuLut_fwd_relpcs_outpcs(
433 icxLuBase *pp,
434 icColorSpaceSignature is, /* Input space, XYZ or Lab */
435 double *out, double *in) {
436 icxLuLut *p = (icxLuLut *)pp;
437
438 /* Convert to the ICC PCS */
439 if (is == icSigLabData && p->natpcs == icSigXYZData) {
440 DBS(("fwd_relpcs_outpcs: Lab in = %s\n", icmPdv(p->inputChan, in)));
441 icmLab2XYZ(&icmD50, out, in);
442 DBS(("fwd_relpcs_outpcs: XYZ = %s\n", icmPdv(p->inputChan, out)));
443 } else if (is == icSigXYZData && p->natpcs == icSigLabData) {
444 DBS(("fwd_relpcs_outpcs: XYZ in = %s\n", icmPdv(p->inputChan, in)));
445 icmXYZ2Lab(&icmD50, out, in);
446 DBS(("fwd_relpcs_outpcs: Lab = %s\n", icmPdv(p->inputChan, out)));
447 } else {
448 DBS(("fwd_relpcs_outpcs: PCS in = %s\n", icmPdv(p->inputChan, in)));
449 icmCpy3(out, in);
450 }
451
452 /* Convert to final PCS (i.e. absolute intent is set that way) */
453 ((icmLuLut *)p->plu)->out_abs((icmLuLut *)p->plu, out, out);
454
455 DBS(("fwd_relpcs_outpcs: abs PCS = %s\n", icmPdv(p->inputChan, out)));
456
457 if (p->outs == icxSigJabData) {
458
459 /* Convert to CAM */
460 p->cam->XYZ_to_cam(p->cam, out, out);
461
462 DBS(("fwd_relpcs_outpcs: Jab = %s\n", icmPdv(p->inputChan, out)));
463 }
464 }
465
466 /* - - - - - - - - - - - - - - - - - - - - - - - - - - */
467 /* Components of inverse lookup, in order */
468
sigfunc(double in,double pp)469 static double sigfunc(double in, double pp) {
470 if (in < 0.5)
471 return pow(2 * in, pp) * 0.5;
472 else
473 return 1.0 - (pow(2 * (1.0 - in), pp) * 0.5);
474 }
475
476 /* Utility function - compute the clip vector direction. */
477 /* return NULL if vector clip isn't used. */
icxClipVector(icxClip * p,double * in,double * cdirv,int safe)478 double *icxClipVector(
479 icxClip *p, /* Clipping setup information */
480 double *in, /* Target point */
481 double *cdirv, /* Returned clip vector */
482 int safe /* Flag - return safe vector */
483 ) {
484 int f;
485 if (p->nearclip != 0)
486 return NULL; /* Doing nearest clipping, not vector */
487
488 if (p->cm == NULL) {
489 /* Default is simple vector clip */
490 for (f = 0; f < p->fdi; f++)
491 cdirv[f] = p->ocent[f] - in[f]; /* Clip towards output gamut center */
492
493 /* Else use CuspMap for more targeted vector */
494 } else {
495 double Cpow = 2.0; /* [2.0] 4.0 for less L change */
496 double Lsig = 2.5; /* [2.5] 1.6 for less L change */
497 double Lpow = 0.5; /* [0.5] 0.8 for less L change */
498 double Cratio = 0.9; /* [0.9] see cusptest.c */
499 double ratio;
500 double ss[3];
501 double inC;
502 double targ[3], cuspLCh[3];
503
504 inC = sqrt(in[1] * in[1] + in[2] * in[2]);
505
506 p->cm->getCusp(p->cm, cuspLCh, in);
507
508
509 //icmLCh2Lab(targ, cuspLCh);
510 //printf("\n~1 in %f %f %f -> cusp %f %f %f\n", in[0], in[1], in[2], targ[0], targ[1], targ[2]);
511 //printf("~1 in %f %f %f -> cuspLCh %f %f %f\n", in[0], in[1], in[2], cuspLCh[0], cuspLCh[1], cuspLCh[2]);
512
513 /* Constrain clip vector to always be inwards pointing */
514 if (cuspLCh[1] > 0.90 * inC) {
515 cuspLCh[1] = 0.90 * inC;
516 //printf("~1 Constrain cusp C: cuspLCh %f %f %f\n", cuspLCh[0], cuspLCh[1], cuspLCh[2]);
517 }
518
519 ss[0] = in[0];
520 ss[1] = in[1];
521 ss[2] = in[2];
522
523 if (ss[0] < p->cm->Lmin[0])
524 ss[0] = p->cm->Lmin[0];
525 if (ss[0] > p->cm->Lmax[0])
526 ss[0] = p->cm->Lmax[0];
527
528 if (safe) {
529 targ[0] = cuspLCh[0]; /* L target is cusp L */
530 targ[1] = 0.0; /* Right on the neutral axis */
531
532 } else {
533 if (ss[0] >= cuspLCh[0]) {
534 ratio = (p->cm->Lmax[0] - ss[0])/(p->cm->Lmax[0] - cuspLCh[0]);
535 targ[0] = p->cm->Lmax[0] - sigfunc(pow(ratio, Lpow), Lsig)
536 * (p->cm->Lmax[0] - cuspLCh[0]);
537 targ[1] = pow(ratio, Cpow) * Cratio * cuspLCh[1];
538 } else {
539 ratio = (ss[0] - p->cm->Lmin[0])/(cuspLCh[0] - p->cm->Lmin[0]);
540 targ[0] = p->cm->Lmin[0] + sigfunc(pow(ratio, Lpow), Lsig)
541 * (cuspLCh[0] - p->cm->Lmin[0]);
542 targ[1] = pow(ratio, Cpow) * Cratio * cuspLCh[1];
543 }
544 }
545 targ[2] = cuspLCh[2]; /* h of in */
546
547 //printf("~1 targLCH %f %f %f\n", targ[0], targ[1], targ[2]);
548 icmLCh2Lab(targ, targ);
549 //printf("~1 targ %f %f %f\n", targ[0], targ[1], targ[2]);
550
551 /* line target point up with white and black point ? */
552 ratio = (ss[0] - p->cm->Lmin[0])/(p->cm->Lmax[0] - p->cm->Lmin[0]);
553 targ[1] += (1.0 - ratio) * p->cm->Lmin[1] + ratio * p->cm->Lmax[1];
554 targ[2] += (1.0 - ratio) * p->cm->Lmin[2] + ratio * p->cm->Lmax[2];
555
556 //printf("~1 in %f %f %f -> targ %f %f %f\n", in[0], in[1], in[2], targ[0], targ[1], targ[2]);
557
558 /* Compute target clip direction */
559 for (f = 0; f < p->fdi; f++)
560 cdirv[f] = (targ[f] - in[f]);
561 }
562
563 return cdirv;
564 }
565
566 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
567 /* Typically inv is used to invert a device->PCS table, */
568 /* so it is doing a PCS->device conversion. */
569 /* This doesn't always have to be the case though. */
570
571 /* PCS override (Effective PCS) to PCS conversion + absolute to relative conversion */
icxLuLut_inv_out_abs(icxLuLut * p,double * out,double * in)572 int icxLuLut_inv_out_abs(icxLuLut *p, double *out, double *in) {
573 int rv = 0;
574
575 DBR(("\nicxLuLut_inv_out_abs got PCS %s\n",icmPdv(p->outputChan, in)));
576
577 if (p->mergeclut == 0) {
578 if (p->outs == icxSigJabData) {
579 p->cam->cam_to_XYZ(p->cam, out, in);
580 DBR(("icxLuLut_inv_out_abs after cam2XYZ %s\n",icmPdv(p->outputChan, out)));
581 /* Hack to prevent CAM02 weirdness being amplified by inv_out_abs() */
582 /* or per channel clipping. */
583 /* Limit -Y to non-stupid values by scaling */
584 if (out[1] < -0.1) {
585 out[0] *= -0.1/out[1];
586 out[2] *= -0.1/out[1];
587 out[1] = -0.1;
588 DBR(("icxLuLut_inv_out_abs after clipping -Y %s\n",icmPdv(p->outputChan, out)));
589 }
590 rv |= ((icmLuLut *)p->plu)->inv_out_abs((icmLuLut *)p->plu, out, out);
591 DBR(("icxLuLut_inv_out_abs after icmLut inv_out_abs %s\n",icmPdv(p->outputChan, out)));
592 } else {
593 rv |= ((icmLuLut *)p->plu)->inv_out_abs((icmLuLut *)p->plu, out, in);
594 DBR(("icxLuLut_inv_out_abs after icmLut inv_out_abs %s\n",icmPdv(p->outputChan, out)));
595 }
596 } else {
597 int i;
598 for (i = 0; i < p->outputChan; i++)
599 out[i] = in[i];
600 }
601 DBR(("icxLuLut_inv_out_abs returning PCS %f %f %f\n",out[0],out[1],out[2]))
602 return rv;
603 }
604
605 /* Do output->output' inverse lookup */
icxLuLut_inv_output(icxLuLut * p,double * out,double * in)606 int icxLuLut_inv_output(icxLuLut *p, double *out, double *in) {
607 int rv = 0;
608 DBR(("icxLuLut_inv_output got PCS %f %f %f\n",in[0],in[1],in[2]))
609 if (p->mergeclut == 0) {
610 #ifdef NEVER
611 rv = ((icmLuLut *)p->plu)->inv_output((icmLuLut *)p->plu, out, in);
612 #else
613 int i,j;
614 int nsoln; /* Number of solutions found */
615 co pp[MAX_INVSOLN]; /* Room for all the solutions found */
616 double cdir;
617
618 for (i = 0; i < p->outputChan; i++) {
619 pp[0].p[0] = p->outputClipc[i];
620 pp[0].v[0] = in[i];
621 cdir = p->outputClipc[i] - in[i]; /* Clip towards output range */
622
623 nsoln = p->outputTable[i]->rev_interp (
624 p->outputTable[i], /* this */
625 RSPL_NEARCLIP, /* Clip to nearest (faster than vector) */
626 MAX_INVSOLN, /* Maximum number of solutions allowed for */
627 NULL, /* No auxiliary input targets */
628 &cdir, /* Clip vector direction and length */
629 pp); /* Input and output values */
630
631 if (nsoln & RSPL_DIDCLIP)
632 rv = 1;
633
634 nsoln &= RSPL_NOSOLNS; /* Get number of solutions */
635
636 if (nsoln == 1) { /* Exactly one solution */
637 j = 0;
638 } else if (nsoln == 0) { /* Zero solutions. This is unexpected. */
639 error("xlut: Unexpected failure to find reverse solution for output table");
640 return 2;
641 } else { /* Multiple solutions */
642 /* Use a simple minded resolution - choose the one closest to the center */
643 double bdist = 1e300;
644 int bsoln = 0;
645 /* Don't expect this - 1D luts are meant to be monotonic */
646 warning("1D lut inversion got %d reverse solutions\n",nsoln);
647 warning("solution 0 = %f\n",pp[0].p[0]);
648 warning("solution 1 = %f\n",pp[1].p[0]);
649 for (j = 0; j < nsoln; j++) {
650 double tt;
651 tt = pp[i].p[0] - p->outputClipc[i];
652 tt *= tt;
653 if (tt < bdist) { /* Better solution */
654 bdist = tt;
655 bsoln = j;
656 }
657 }
658 j = bsoln;
659 }
660 out[i] = pp[j].p[0];
661 }
662
663 #endif /* NEVER */
664 } else {
665 int i;
666 for (i = 0; i < p->outputChan; i++)
667 out[i] = in[i];
668 }
669 DBR(("icxLuLut_inv_output returning PCS' %f %f %f\n",out[0],out[1],out[2]))
670 return rv;
671 }
672
673 /* Ink limit+gamut limit calculation function for xLuLut. */
674 /* Returns < 0.0 if input value is within limits, */
675 /* > 0.0 if outside limits. Value is proportinal to distance to limits. */
676 /* We implement device gamut check to improve utility outside rspl, */
677 /* in optimisation routines. */
678 /* The limits are assumed to be post calibrated device values (ie. real */
679 /* final device values) */
icxLimit(icxLuLut * p,double * in)680 static double icxLimit(
681 icxLuLut *p,
682 double *in
683 ) {
684 double cin[MAX_CHAN]; /* Calibrated input values */
685 double tlim, klim;
686 double ovr, val;
687 int e;
688
689 if (p->pp->cal != NULL) { /* We have device calibration information */
690 p->pp->cal->interp(p->pp->cal, cin, in);
691 } else {
692 for (e = 0; e < p->inputChan; e++)
693 cin[e] = in[e];
694 }
695
696 if ((tlim = p->ink.tlimit) < 0.0)
697 tlim = (double)p->inputChan; /* Default is no limit */
698
699 if ((klim = p->ink.klimit) < 0.0)
700 klim = 1.0;
701
702 /* Compute amount outside total limit */
703 { /* No calibration */
704 double sum;
705 for (sum = 0.0, e = 0; e < p->inputChan; e++)
706 sum += cin[e];
707 val = sum - tlim;
708 }
709
710 /* Compute amount outside black limit */
711 if (p->ink.klimit >= 0.0) {
712 double kval = 0.0;
713 switch(p->natis) {
714 case icSigCmykData:
715 kval = cin[3] - klim;
716 break;
717 default:
718 /* NOTE !!! p->kch isn't being initialized !!! */
719 if (p->kch >= 0) {
720 kval = cin[p->kch] - klim;
721 } else {
722 error("xlut: Unknown colorspace when black limit specified");
723 }
724 }
725 if (kval > val)
726 val = kval;
727 }
728 /* Compute amount outside device value limits 0.0 - 1.0 */
729 for (ovr = -1.0, e = 0; e < p->inputChan; e++) {
730 if (in[e] < 0.0) { /* ! Not cin[] */
731 if (-in[e] > ovr)
732 ovr = -in[e];
733 } else if (in[e] > 1.0) {
734 if ((in[e] - 1.0) > ovr)
735 ovr = in[e] - 1.0;
736 }
737 }
738 if (ovr > val)
739 val = ovr;
740
741 return val;
742 }
743
744 /* Same as above, but works with input' values */
745 /* (If an ink limit is being used we assume that the */
746 /* input space is not PCS, hence inv_in_abs() is doing nothing) */
icxLimitD(icxLuLut * p,double * ind)747 static double icxLimitD(
748 icxLuLut *p,
749 double *ind
750 ) {
751 double in[MAX_CHAN];
752 co tc;
753 int e;
754
755 /* Convert input' to input through revinput Luts (for speed) */
756 for (e = 0; e < p->inputChan; e++) {
757 tc.p[0] = ind[e];
758 p->revinputTable[e]->interp(p->revinputTable[e], &tc);
759 in[e] = tc.v[0];
760 }
761
762 return icxLimit(p, in);
763 }
764
765 /* Ink limit+gamut limit clipping function for xLuLut (CMYK). */
766 /* Return nz if there was clipping */
icxDoLimit(icxLuLut * p,double * out,double * in)767 static int icxDoLimit(
768 icxLuLut *p,
769 double *out,
770 double *in
771 ) {
772 double tlim, klim = -1.0;
773 double sum;
774 int clip = 0, e;
775 int kch = -1;
776
777 for (e = 0; e < p->inputChan; e++)
778 out[e] = in[e];
779
780 if ((tlim = p->ink.tlimit) < 0.0)
781 tlim = (double)p->inputChan;
782
783 if ((klim = p->ink.klimit) < 0.0)
784 klim = 1.0;
785
786 /* Clip black */
787 if (p->natis == icSigCmykData)
788 kch = 3;
789 else
790 kch = p->kch;
791
792 if (kch >= 0) {
793 if (out[p->kch] > klim) {
794 out[p->kch] = klim;
795 clip = 1;
796 }
797 }
798
799 /* Compute amount outside total limit */
800 for (sum = 0.0, e = 0; e < p->inputChan; e++)
801 sum += out[e];
802
803 if (sum > tlim) {
804 clip = 1;
805 sum /= (double)p->inputChan;
806 for (e = 0; e < p->inputChan; e++)
807 out[e] -= sum;
808 }
809 return clip;
810 }
811
812 #ifdef NEVER
813 #undef DBK
814 #define DBK(xxx) printf xxx ;
815 #else
816 #undef DBK
817 #define DBK(xxx)
818 #endif
819
820 /* helper function that creates our standard K locus curve value, */
821 /* given the curve parameters, and the normalised L 0.0 - 1.0 value. */
822 /* No filtering version. */
823 /* !!! Should add K limit in here so that smoothing takes it into account !!! */
icxKcurveNF(double L,icxInkCurve * c)824 static double icxKcurveNF(double L, icxInkCurve *c) {
825 double Kstpo, Kenpo, Kstle, Kenle;
826 double rv;
827
828 DBK(("icxKcurve got L = %f, smth %f skew %f, Parms %f %f %f %f %f\n",L, c->Ksmth, c->Kskew, c->Kstle, c->Kstpo, c->Kenpo, c->Kenle, c->Kshap));
829
830 /* Invert sense of L, so that 0.0 = white, 1.0 = black */
831 L = 1.0 - L;
832
833 /* Clip L, just in case */
834 if (L < 0.0) {
835 L = 0.0;
836 } else if (L > 1.0) {
837 L = 1.0;
838 }
839 DBK(("Clipped inverted L = %f\n",L));
840
841 /* Make sure breakpoints are ordered */
842 if (c->Kstpo < c->Kenpo) {
843 Kstle = c->Kstle;
844 Kstpo = c->Kstpo;
845 Kenpo = c->Kenpo;
846 Kenle = c->Kenle;
847 } else { /* They're swapped */
848 Kstle = c->Kenle;
849 Kstpo = c->Kenpo;
850 Kenpo = c->Kstpo;
851 Kenle = c->Kstle;
852 }
853
854 if (L <= Kstpo) {
855 /* We are at white level */
856 rv = Kstle;
857 DBK(("At white level %f\n",rv));
858 } else if (L >= Kenpo) {
859 /* We are at black level */
860 rv = Kenle;
861 DBK(("At black level %f\n",rv));
862 } else {
863 double Lp, g;
864 /* We must be on the curve from start to end levels */
865
866 Lp = (L - Kstpo)/(Kenpo - Kstpo);
867
868 DBK(("Curve position %f\n",Lp));
869
870 Lp = pow(Lp, c->Kskew);
871
872 DBK(("Skewed curve position %f\n",Lp));
873
874 g = c->Kshap/2.0;
875 DBK(("Curve bf %f, g %g\n",Lp,g));
876
877 /* A value of 0.5 will be tranlated to g */
878 Lp = Lp/((1.0/g - 2.0) * (1.0 - Lp) + 1.0);
879
880 DBK(("Skewed shaped %f\n",Lp));
881
882 Lp = pow(Lp, 1.0/c->Kskew);
883
884 DBK(("Shaped %f\n",Lp));
885
886 /* Transition between start end end levels */
887 rv = Lp * (Kenle - Kstle) + Kstle;
888
889 DBK(("Scaled to start and end levele %f\n",rv));
890 }
891
892 DBK(("Returning %f\n",rv));
893 return rv;
894 }
895
896 #ifdef DBK
897 #undef DBK
898 #define DBK(xxx)
899 #endif
900
901
902 #ifdef NEVER
903 #undef DBK
904 #define DBK(xxx) printf xxx ;
905 #else
906 #undef DBK
907 #define DBK(xxx)
908 #endif
909
910 /* Same as above, but implement transition filters across inflection points. */
911 /* (The convolultion filter previously used could be */
912 /* re-instituted if something was done about compressing */
913 /* the filter at the boundaries so that the levels are met.) */
icxKcurve(double L,icxInkCurve * c)914 static double icxKcurve(double L, icxInkCurve *c) {
915
916 #ifdef DISABLE_KCURVE_FILTER
917 return icxKcurveNF(L, c);
918
919 #else /* !DISABLE_KCURVE_FILTER */
920
921 double Kstpo, Kenpo, Kstle, Kenle, Ksmth;
922 double rv;
923
924 DBK(("icxKcurve got L = %f, smth %f skew %f, Parms %f %f %f %f %f\n",L, c->Ksmth, c->Kskew, c->Kstle, c->Kstpo, c->Kenpo, c->Kenle, c->Kshap));
925
926 /* Invert sense of L, so that 0.0 = white, 1.0 = black */
927 L = 1.0 - L;
928
929 /* Clip L, just in case */
930 if (L < 0.0) {
931 L = 0.0;
932 } else if (L > 1.0) {
933 L = 1.0;
934 }
935 DBK(("Clipped inverted L = %f\n",L));
936
937 /* Make sure breakpoints are ordered */
938 if (c->Kstpo < c->Kenpo) {
939 Kstle = c->Kstle;
940 Kstpo = c->Kstpo;
941 Kenpo = c->Kenpo;
942 Kenle = c->Kenle;
943 } else { /* They're swapped */
944 Kstle = c->Kenle;
945 Kstpo = c->Kenpo;
946 Kenpo = c->Kstpo;
947 Kenle = c->Kstle;
948 }
949 Ksmth = c->Ksmth;
950
951 /* Curve value at point */
952 rv = icxKcurveNF(1.0 - L, c);
953
954 DBK(("Raw output at iL = %f, rv\n",L));
955
956 /* Create filtered value */
957 {
958 double wbs, wbe; /* White transitioin start, end */
959 double wbl, wfv; /* White blend factor, value at filter band */
960
961 double mt; /* Middle of the two transitions */
962
963 double bbs, bbe; /* Black transitioin start, end */
964 double bbl, bfv; /* Black blend factor, value at filter band */
965
966 wbs = Kstpo - Ksmth;
967 wbe = Kstpo + Ksmth;
968
969 bbs = Kenpo - 1.0 * Ksmth;
970 bbe = Kenpo + 1.0 * Ksmth;
971
972 mt = 0.5 * (wbe + bbs);
973
974 /* Make sure that the whit & black transition regions */
975 /* don't go out of bounts or overlap */
976 if (wbs < 0.0) {
977 wbe += wbs;
978 wbs = 0.0;
979 }
980 if (bbe > 1.0) {
981 bbs += (bbe - 1.0);
982 bbe = 1.0;
983 }
984
985 if (wbe > mt) {
986 wbs += (wbe - mt);
987 wbe = mt;
988 }
989
990 if (bbs < mt) {
991 bbe += (mt - bbs);
992 bbs = mt;
993 }
994
995 DBK(("Transition windows %f - %f, %f - %f\n",wbs, wbe, bbw, bbe));
996 if (wbs < wbe) {
997 wbl = (L - wbe)/(wbs - wbe);
998
999 if (wbl > 0.0 && wbl < 1.0) {
1000 wfv = icxKcurveNF(1.0 - wbe, c);
1001 DBK(("iL = %f, wbl = %f, wfv = %f\n",L,Kstpo,wbl,wfv));
1002
1003 wbl = 1.0 - pow(1.0 - wbl, 2.0);
1004 rv = wbl * Kstle + (1.0 - wbl) * wfv;
1005 }
1006 }
1007 if (bbs < bbe) {
1008 bbl = (L - bbe)/(bbs - bbe);
1009
1010 if (bbl > 0.0 && bbl < 1.0) {
1011 bfv = icxKcurveNF(1.0 - bbs, c);
1012 DBK(("iL = %f, bbl = %f, bfv = %f\n",L,Kstpo,bbl,bfv));
1013
1014 bbl = pow(bbl, 2.0);
1015 rv = bbl * bfv + (1.0 - bbl) * Kenle;
1016 }
1017 }
1018 }
1019
1020 /* To be safe */
1021 if (rv < 0.0)
1022 rv = 0.0;
1023 else if (rv > 1.0)
1024 rv = 1.0;
1025
1026 DBK(("Returning %f\n",rv));
1027 return rv;
1028 #endif /* !DISABLE_KCURVE_FILTER */
1029 }
1030
1031 #ifdef DBK
1032 #undef DBK
1033 #define DBK(xxx)
1034 #endif
1035
1036 /* Do output'->input' lookup with aux details. */
1037 /* Note that out[] will be used as the inking value if icxKrule is */
1038 /* icxKvalue, icxKlocus, icxKl5l or icxKl5lk, and that the auxiliar values, PCS ranges */
1039 /* and icxKrule value will all be evaluated in output->input space (not ' space). */
1040 /* Note that the ink limit will be computed after converting input' to input, auxt */
1041 /* will override the inking rule, and auxr[] reflects the available auxiliary range */
1042 /* that the locus was to choose from, and auxv[] was the actual auxiliary used. */
1043 /* Returns clip status. */
icxLuLut_inv_clut_aux(icxLuLut * p,double * out,double * auxv,double * auxr,double * auxt,double * clipd,double * in)1044 int icxLuLut_inv_clut_aux(
1045 icxLuLut *p,
1046 double *out, /* Function return values, plus aux value or locus target input if auxt == NULL */
1047 double *auxv, /* If not NULL, return aux value used (packed) */
1048 double *auxr, /* If not NULL, return aux locus range (packed, 2 at a time) */
1049 double *auxt, /* If not NULL, specify the aux target for this lookup (override ink) */
1050 double *clipd, /* If not NULL, return DE to gamut on clipi, 0 for not clip */
1051 double *in /* Function input values to invert (== clut output' values) */
1052 ) {
1053 co pp[MAX_INVSOLN]; /* Room for all the solutions found */
1054 co upp; /* pp[0] value sent to rev_interp() for replay. */
1055 int nsoln; /* Number of solutions found */
1056 double *cdir, cdirv[MXDO]; /* Clip vector direction and length/LCh weighting */
1057 int e,f,i;
1058 int fdi = p->clutTable->fdi;
1059 int flags = 0; /* reverse interp flags */
1060 int xflags = 0; /* extra clip/noclip flags */
1061 int uflags = 0; /* flags value sent to rev_interp() for replay */
1062 double tin[MXDO]; /* PCS value to be inverted */
1063 double cdist = 0.0; /* clip DE */
1064 int crv = 0; /* Return value - set to 1 if clipped */
1065
1066 if (p->nearclip != 0)
1067 flags |= RSPL_NEARCLIP; /* Use nearest clipping rather than clip vector */
1068
1069 DBR(("inv_clut_aux input is %f %f %f\n",in[0], in[1], in[2]))
1070
1071 if (auxr != NULL) { /* Set a default locus range */
1072 int ee = 0;
1073 for (e = 0; e < p->clutTable->di; e++) {
1074 if (p->auxm[e] != 0) {
1075 auxr[ee++] = 1e60;
1076 auxr[ee++] = -1e60;
1077 }
1078 }
1079 }
1080
1081 /* Setup for reverse lookup */
1082 for (f = 0; f < fdi; f++)
1083 upp.v[f] = pp[0].v[f] = in[f]; /* Target value */
1084
1085 /* Compute clip vector, if any */
1086 cdir = icxClipVector(&p->clip, in, cdirv, 0);
1087
1088 if (p->clutTable->di > fdi) { /* ie. CMYK->Lab, there will be ambiguity */
1089 double min[MXDI], max[MXDI]; /* Auxiliary locus range */
1090
1091 #ifdef REPORT_LOCUS_SEGMENTS /* Examine how many segments there are */
1092 { /* ie. CMYK->Lab, there will be ambiguity */
1093 double smin[10][MXRI], smax[10][MXRI]; /* Auxiliary locus segment ranges */
1094 double min[MXRI], max[MXRI]; /* Auxiliary min and max locus range */
1095
1096 nsoln = p->clutTable->rev_locus_segs(
1097 p->clutTable, /* rspl object */
1098 p->auxm, /* Auxiliary mask */
1099 pp, /* Input target and output solutions */
1100 10, /* Maximum number of solutions to return */
1101 smin, smax); /* Returned locus of valid auxiliary values */
1102
1103 if (nsoln != 0) {
1104 /* Convert the locuses from input' -> input space */
1105 /* and get overall min/max locus range */
1106 for (e = 0; e < p->clutTable->di; e++) {
1107 co tc;
1108 /* (Is speed more important than precision ?) */
1109 if (p->auxm[e] != 0) {
1110 for (i = 0; i < nsoln; i++) {
1111 tc.p[0] = smin[i][e];
1112 p->revinputTable[e]->interp(p->revinputTable[e], &tc);
1113 smin[i][e] = tc.v[0];
1114 tc.p[0] = smax[i][e];
1115 p->revinputTable[e]->interp(p->revinputTable[e], &tc);
1116 smax[i][e] = tc.v[0];
1117 printf(" Locus seg %d:[%d] %f -> %f\n",i, e, smin[i][e], smax[i][e]);
1118 }
1119 }
1120 }
1121 }
1122 }
1123 #endif /* REPORT_LOCUS_SEGMENTS */
1124
1125 /* Compute auxiliary locus on the fly. This is in dev' == input' space. */
1126 nsoln = p->clutTable->rev_locus(
1127 p->clutTable, /* rspl object */
1128 p->auxm, /* Auxiliary mask */
1129 pp, /* Input target and output solutions */
1130 min, max); /* Returned locus of valid auxiliary values */
1131
1132 if (nsoln == 0) {
1133 xflags |= RSPL_WILLCLIP; /* No valid locus, so we expect to have to clip */
1134 #ifdef DEBUG_RLUT
1135 printf("inv_clut_aux: no valid locus, expect clip\n");
1136 #endif
1137 /* Make sure that the auxiliar value is initialized, */
1138 /* even though it won't have any effect. */
1139 for (e = 0; e < p->clutTable->di; e++) {
1140 if (p->auxm[e] != 0) {
1141 upp.p[e] = pp[0].p[e] = 0.5;
1142 }
1143 }
1144
1145 } else { /* Got a valid locus */
1146
1147 /* Convert the locuses from input' -> input space */
1148 for (e = 0; e < p->clutTable->di; e++) {
1149 co tc;
1150 /* (Is speed more important than precision ?) */
1151 if (p->auxm[e] != 0) {
1152 tc.p[0] = min[e];
1153 p->revinputTable[e]->interp(p->revinputTable[e], &tc);
1154 min[e] = tc.v[0];
1155 tc.p[0] = max[e];
1156 p->revinputTable[e]->interp(p->revinputTable[e], &tc);
1157 max[e] = tc.v[0];
1158 }
1159 }
1160
1161 if (auxr != NULL) { /* Report the locus range */
1162 int ee = 0;
1163 for (e = 0; e < p->clutTable->di; e++) {
1164 if (p->auxm[e] != 0) {
1165 auxr[ee++] = min[e];
1166 auxr[ee++] = max[e];
1167 }
1168 }
1169 }
1170
1171 if (auxt != NULL) { /* overiding auxiliary target */
1172 int ee = 0;
1173 for (e = 0; e < p->clutTable->di; e++) {
1174 if (p->auxm[e] != 0) {
1175 double iv = auxt[ee++];
1176 if (iv < min[e])
1177 iv = min[e];
1178 else if (iv > max[e])
1179 iv = max[e];
1180 upp.p[e] = pp[0].p[e] = iv;
1181 }
1182 }
1183 DBR(("inv_clut_aux: aux %f from auxt[] %f\n",pp[0].p[3],auxt[0]))
1184 } else if (p->ink.k_rule == icxKvalue) {
1185 /* Implement the auxiliary inking rule */
1186 /* Target auxiliary values are provided in out[] K value */
1187 for (e = 0; e < p->clutTable->di; e++) {
1188 if (p->auxm[e] != 0) {
1189 double iv = out[e]; /* out[] holds aux target value */
1190 if (iv < min[e])
1191 iv = min[e];
1192 else if (iv > max[e])
1193 iv = max[e];
1194 upp.p[e] = pp[0].p[e] = iv;
1195 }
1196 }
1197 DBR(("inv_clut_aux: aux %f from out[0] K target %f min %f max %f\n",pp[0].p[3],out[3],min[3],max[3]))
1198 } else if (p->ink.k_rule == icxKlocus) {
1199 /* Set target auxliary input values from values in out[] and locus */
1200 for (e = 0; e < p->clutTable->di; e++) {
1201 if (p->auxm[e] != 0) {
1202 double ii, iv;
1203 ii = out[e]; /* Input ink locus */
1204 iv = min[e] + ii * (max[e] - min[e]); /* Output ink from locus */
1205 if (iv < min[e])
1206 iv = min[e];
1207 else if (iv > max[e])
1208 iv = max[e];
1209 upp.p[e] = pp[0].p[e] = iv;
1210 }
1211 }
1212 DBR(("inv_clut_aux: aux %f from out[0] locus %f min %f max %f\n",pp[0].p[3],out[3],min[3],max[3]))
1213 } else { /* p->ink.k_rule == icxKluma5 || icxKluma5k || icxKl5l || icxKl5lk */
1214 /* Auxiliaries are driven by a rule and the output values */
1215 double rv, L;
1216
1217 /* If we've got a mergeclut, then the PCS' is the same as the */
1218 /* effective PCS, and we need to convert to native PCS */
1219 if (p->mergeclut) {
1220 p->mergeclut = 0; /* Hack to be able to use inv_out_abs() */
1221 icxLuLut_inv_out_abs(p, tin, in);
1222 p->mergeclut = 1;
1223
1224 } else {
1225 /* Convert native PCS' to native PCS values */
1226 p->output(p, tin, in);
1227 }
1228
1229 /* Figure out Luminance number */
1230 if (p->natos == icSigXYZData) {
1231 icmXYZ2Lab(&icmD50, tin, tin);
1232 } else if (p->natos != icSigLabData) { /* Hmm. that's unexpected */
1233 error("Assert: xlut K locus, unexpected native pcs of 0x%x\n",p->natos);
1234 }
1235 L = 0.01 * tin[0];
1236 DBR(("inv_clut_aux: aux from Luminance, raw L = %f\n",L));
1237
1238 /* Normalise L to its possible range from min to max */
1239 L = (L - p->Lmin)/(p->Lmax - p->Lmin);
1240 DBR(("inv_clut_aux: Normalize L = %f\n",L));
1241
1242 /* Convert L to curve value */
1243 rv = icxKcurve(L, &p->ink.c);
1244 DBR(("inv_clut_aux: icxKurve lookup returns = %f\n",rv));
1245
1246 if (p->ink.k_rule == icxKluma5) { /* Curve is locus value */
1247
1248 /* Set target black as K fraction within locus */
1249
1250 for (e = 0; e < p->clutTable->di; e++) {
1251 if (p->auxm[e] != 0) {
1252 upp.p[e] = pp[0].p[e] = min[e] + rv * (max[e] - min[e]);
1253 }
1254 }
1255 DBR(("inv_clut_aux: aux %f from locus %f min %f max %f\n",pp[0].p[3],rv,min[3],max[3]))
1256
1257 } else if (p->ink.k_rule == icxKluma5k) { /* Curve is K value */
1258
1259 for (e = 0; e < p->clutTable->di; e++) {
1260 if (p->auxm[e] != 0) {
1261 double iv = rv;
1262 if (iv < min[e]) /* Clip to locus */
1263 iv = min[e];
1264 else if (iv > max[e])
1265 iv = max[e];
1266 upp.p[e] = pp[0].p[e] = iv;
1267 }
1268 }
1269 DBR(("inv_clut_aux: aux %f from out[0] K target %f min %f max %f\n",pp[0].p[3],rv,min[3],max[3]))
1270
1271 } else { /* icxKl5l || icxKl5lk */
1272 /* Create second curve, and use input locus to */
1273 /* blend between */
1274
1275 double rv2; /* Upper limit */
1276
1277 /* Convert L to max curve value */
1278 rv2 = icxKcurve(L, &p->ink.x);
1279
1280 if (rv2 < rv) { /* Ooops - better swap. */
1281 double tt;
1282 tt = rv;
1283 rv = rv2;
1284 rv2 = tt;
1285 }
1286
1287 for (e = 0; e < p->clutTable->di; e++) {
1288 if (p->auxm[e] != 0) {
1289 if (p->ink.k_rule == icxKl5l) {
1290 double ii;
1291 ii = out[e]; /* Input K locus */
1292 if (ii < 0.0)
1293 ii = 0.0;
1294 else if (ii > 1.0)
1295 ii = 1.0;
1296 ii = (1.0 - ii) * rv + ii * rv2;/* Blend between locus rule curves */
1297 /* Out ink from output locus */
1298 upp.p[e] = pp[0].p[e] = min[e] + ii * (max[e] - min[e]);
1299 } else {
1300 double iv;
1301 iv = out[e]; /* Input K level */
1302 if (iv < rv) /* Constrain to curves */
1303 iv = rv;
1304 else if (iv > rv2)
1305 iv = rv2;
1306 upp.p[e] = pp[0].p[e] = iv;
1307 }
1308 }
1309 }
1310 DBR(("inv_clut_aux: aux %f from 2 curves\n",pp[0].p[3]))
1311 }
1312 }
1313
1314 /* Convert to input/dev aux target to input'/dev' space for rspl inversion */
1315 for (e = 0; e < p->clutTable->di; e++) {
1316 double tv, bv = 0.0, bd = 1e6;
1317 co tc;
1318 if (p->auxm[e] != 0) {
1319 tv = pp[0].p[e];
1320 /* Clip to overall locus range (belt and braces) */
1321 if (tv < min[e])
1322 tv = min[e];
1323 if (tv > max[e])
1324 tv = max[e];
1325 tc.p[0] = tv;
1326 p->inputTable[e]->interp(p->inputTable[e], &tc);
1327 upp.p[e] = pp[0].p[e] = tc.v[0];
1328 }
1329 }
1330
1331 xflags |= RSPL_EXACTAUX; /* Since we confine aux to locus */
1332
1333 #ifdef DEBUG_RLUT
1334 printf("inv_clut_aux computed aux values ");
1335 for (e = 0; e < p->clutTable->di; e++) {
1336 if (p->auxm[e] != 0)
1337 printf("%d: %f ",e,pp[0].p[e]);
1338 }
1339 printf("\n");
1340 #endif /* DEBUG_RLUT */
1341 }
1342
1343 if (clipd != NULL) { /* Copy pp.v[] to compute clip DE */
1344 for (f = 0; f < fdi; f++)
1345 tin[f] = pp[0].v[f];
1346 }
1347
1348 uflags = RSPL_MAXAUX | flags | xflags; /* Combine all the flags */
1349
1350 /* Find reverse solution with target auxiliaries */
1351 /* We choose the closest aux at or above the target */
1352 /* to try and avoid glitches near black due to */
1353 /* possible forked black locuses. */
1354 nsoln = p->clutTable->rev_interp(
1355 p->clutTable, /* rspl object */
1356 uflags,
1357 MAX_INVSOLN, /* Maxumum solutions to return */
1358 p->auxm, /* Auxiliary input chanel mask */
1359 cdir, /* Clip vector direction/LCh weighting */
1360 pp); /* Input target and output solutions */
1361 /* returned solutions in pp[0..retval-1].p[] */
1362
1363 } else {
1364 DBR(("inv_clut_aux needs no aux value\n"))
1365
1366 if (clipd != NULL) { /* Copy pp.v[] to compute clip DE */
1367 for (f = 0; f < fdi; f++)
1368 tin[f] = pp[0].v[f];
1369 }
1370
1371 uflags = flags; /* No extra flags */
1372
1373 /* Color spaces don't need auxiliaries to choose from solution locus */
1374 nsoln = p->clutTable->rev_interp(
1375 p->clutTable, /* rspl object */
1376 uflags,
1377 MAX_INVSOLN, /* Maxumum solutions to return */
1378 NULL, /* No auxiliary input targets */
1379 cdir, /* Clip vector direction/LCh weighting */
1380 pp); /* Input target and output solutions */
1381 /* returned solutions in pp[0..retval-1].p[] */
1382 }
1383 if (nsoln & RSPL_DIDCLIP)
1384 crv = 1; /* Clipped on PCS inverse lookup */
1385
1386 if (crv && clipd != NULL) {
1387
1388
1389 /* Compute the clip DE */
1390 for (cdist = 0.0, f = 0; f < fdi; f++) {
1391 double tt;
1392 tt = pp[0].v[f] - tin[f];
1393 cdist += tt * tt;
1394 }
1395 cdist = sqrt(cdist);
1396 DBR(("Targ PCS %f %f %f Clipped %f %f %f cdist %f\n",tin[0],tin[1],tin[2],pp[0].v[0],pp[0].v[1],pp[0].v[2],cdist))
1397 }
1398
1399 nsoln &= RSPL_NOSOLNS; /* Get number of solutions */
1400
1401 DBR(("inv_clut_aux got %d rev_interp solutions, clipflag = %d\n",nsoln,crv))
1402
1403 /* If we clipped and we should clip in CAM Jab space, compute reverse */
1404 /* clip solution using our additional CAM space rspl. */
1405 /* Note that we don't support vector clip in CAM space at the moment */
1406 /* - icxLuLut_init_clut_camclip() doesn't setup CAM rspl for vector clip. */
1407 if (crv != 0 && p->camclip && p->nearclip) {
1408 co cpp; /* Alternate CAM space solution */
1409 double bf; /* Blend factor */
1410
1411 DBR(("inv_clut_aux got clip, compute CAM clip\n"))
1412
1413 if (nsoln != 1) { /* This would be unexpected */
1414 error("Unexpected failure to return 1 solution on clip for input to output table");
1415 }
1416
1417 if (p->cclutTable == NULL) { /* we haven't created this yet, so do so */
1418 if (icxLuLut_init_clut_camclip(p))
1419 error("Creating CAM rspl for camclip failed");
1420 }
1421
1422 /* Setup for reverse lookup */
1423 DBR(("inv_clut_aux cam clip PCS in %f %f %f\n",in[0],in[1],in[2]))
1424
1425 /* Convert from PCS' to (XYZ) PCS */
1426 ((icmLuLut *)p->absxyzlu)->output((icmLuLut *)p->absxyzlu, tin, in);
1427 DBR(("inv_clut_aux cam clip PCS' -> PCS %f %f %f\n",tin[0],tin[1],tin[2]))
1428
1429 ((icmLuLut *)p->absxyzlu)->out_abs((icmLuLut *)p->absxyzlu, tin, tin);
1430 DBR(("inv_clut_aux cam clip abs XYZ PCS %f %f %f\n",tin[0],tin[1],tin[2]))
1431
1432 p->cam->XYZ_to_cam(p->cam, tin, tin);
1433 DBR(("inv_clut_aux cam clip PCS after XYZtoCAM %f %f %f\n",tin[0],tin[1],tin[2]))
1434
1435 for (f = 0; f < fdi; f++) /* Transfer CAM targ */
1436 cpp.v[f] = tin[f];
1437
1438 /* Make sure that the auxiliar value is initialized, */
1439 /* even though it shouldn't have any effect, since should clipp. */
1440 for (e = 0; e < p->clutTable->di; e++) {
1441 if (p->auxm[e] != 0) {
1442 cpp.p[e] = 0.5;
1443 }
1444 }
1445 if (p->clutTable->di > fdi) { /* ie. CMYK->Lab, there will be ambiguity */
1446
1447 nsoln = p->cclutTable->rev_interp(
1448 p->cclutTable, /* rspl object */
1449 flags | xflags | RSPL_WILLCLIP, /* Combine all the flags + clip ?? */
1450 1, /* Maximum solutions to return */
1451 p->auxm, /* Auxiliary input chanel mask */
1452 cdir, /* Clip vector direction/ LCh weighting */
1453 &cpp); /* Input target and output solutions */
1454
1455 } else {
1456 nsoln = p->cclutTable->rev_interp(
1457 p->cclutTable, /* rspl object */
1458 flags | RSPL_WILLCLIP, /* Because we know it will clip ?? */
1459 1, /* Maximum solutions to return */
1460 NULL, /* No auxiliary input targets */
1461 cdir, /* Clip vector direction/LCh weighting */
1462 &cpp); /* Input target and output solutions */
1463 }
1464
1465 nsoln &= RSPL_NOSOLNS; /* Get number of solutions */
1466
1467 if (nsoln != 1) { /* This would be unexpected */
1468 error("Unexpected failure to return 1 solution on CAM clip for input to output table");
1469 }
1470
1471 /* Compute the CAM clip distances */
1472 for (cdist = 0.0, f = 0; f < fdi; f++) {
1473 double tt;
1474 tt = cpp.v[f] - tin[f];
1475 cdist += tt * tt;
1476 }
1477 cdist = sqrt(cdist);
1478 DBR(("Targ CAM %f %f %f Clipped %f %f %f cdist %f\n",tin[0],tin[1],tin[2],cpp.v[0],cpp.v[1],cpp.v[2],cdist))
1479
1480 /* Use magic number to set blend distance, and compute a blend factor. */
1481 /* Blend over 1 delta E, otherwise the Lab & CAM02 divergence can result */
1482 /* in reversals. */
1483 bf = cdist/CAMCLIPTRANS; /* 0.0 for PCS result, 1.0 for CAM result */
1484 if (bf > 1.0)
1485 bf = 1.0;
1486 #ifdef USECAMCLIPSPLINE
1487 bf = bf * bf * (3.0 - 2.0 * bf); /* Convert to spline blend */
1488 #endif
1489 DBR(("cdist %f, spline blend %f\n",cdist,bf))
1490
1491 /* Blend between solution values for PCS and CAM clip result. */
1492 /* We're hoping that the solutions are close, and expect them to be */
1493 /* that way when we're close to the gamut surface (since the cell */
1494 /* vertex values should be exact, irrespective of which interpolation */
1495 /* space they're in), but weird things could happen away from the surface. */
1496 /* Weird things can happen anyway with "clip to nearest", since this is not */
1497 /* guaranteed to be a smooth function, depending on the gamut surface */
1498 /* geometry, without taking some precaution such as clipping to a */
1499 /* convex hull "wrapper" to create a clip vector, which we're not */
1500 /* current doing. */
1501 DBR(("Clip blend between device:\n"))
1502 DBR(("Lab: %s & ",icmPdv(p->clutTable->di, pp[0].p)))
1503 DBR(("Jab: %s ",icmPdv(p->clutTable->di, cpp.p)))
1504
1505 for (e = 0; e < p->clutTable->di; e++) {
1506 out[e] = (1.0 - bf) * pp[0].p[e] + bf * cpp.p[e];
1507 }
1508 DBR((" = %s\n",icmPdv(p->clutTable->di, out)))
1509
1510 /* Not CAM clip case */
1511 } else {
1512
1513 /* If vector clip, replay with simple vector */
1514 if (nsoln == 0 && p->nearclip == 0) {
1515 //printf("~1 !!!!!!!!!!!!!! Vector clip failed - trying safe replay !!!!!!!!!!!!!!!1\n");
1516
1517 // Reset pp[0] target values and auiliaries
1518 for (e = 0; e < p->clutTable->di; e++)
1519 pp[0].p[e] = upp.p[e];
1520 for (f = 0; f < fdi; f++)
1521 pp[0].v[f] = upp.v[f];
1522
1523 /* Compute safe clip vector */
1524 cdir = icxClipVector(&p->clip, in, cdirv, 1);
1525
1526 /* Replay the lookup using safer vector */
1527 nsoln = p->clutTable->rev_interp(
1528 p->clutTable, /* rspl object */
1529 uflags,
1530 MAX_INVSOLN, /* Maxumum solutions to return */
1531 NULL, /* No auxiliary input targets */
1532 cdir, /* Clip vector direction and length */
1533 pp); /* Input target and output solutions */
1534 /* returned solutions in pp[0..retval-1].p[] */
1535 nsoln &= RSPL_NOSOLNS; /* Get number of solutions */
1536
1537 if (nsoln == 0) { /* Hmm */
1538 //printf("~1 !!!!!!!!!!!!!! Vector clip failed again - using nn replay !!!!!!!!!!!!!!!1\n");
1539 // Reset pp[0] target values and auiliaries
1540 for (e = 0; e < p->clutTable->di; e++)
1541 pp[0].p[e] = upp.p[e];
1542 for (f = 0; f < fdi; f++)
1543 pp[0].v[f] = upp.v[f];
1544
1545 /* Replay the lookup as a nearclip, to guarantee a solution */
1546 nsoln = p->clutTable->rev_interp(
1547 p->clutTable, /* rspl object */
1548 RSPL_NEARCLIP | RSPL_NONNSETUP,
1549 MAX_INVSOLN, /* Maxumum solutions to return */
1550 NULL, /* No auxiliary input targets */
1551 NULL, /* No LCh weighting */
1552 pp); /* Input target and output solutions */
1553 /* returned solutions in pp[0..retval-1].p[] */
1554 nsoln &= RSPL_NOSOLNS; /* Get number of solutions */
1555 }
1556 }
1557
1558 if (nsoln == 1) { /* Exactly one solution */
1559 i = 0;
1560 } else if (nsoln == 0) { /* Zero solutions. This is unexpected. */
1561 double in_v[MXDO];
1562 p->output(p, in_v, pp[0].v); /* Get ICC inverse input values */
1563 p->out_abs(p, in_v, in_v);
1564 if (p->nearclip == 0)
1565 a1logd(g_log,0,"Clip dst %f %f %f\n",pp[0].v[0]+cdir[0], pp[0].v[1]+cdir[1], pp[0].v[2]+cdir[2]);
1566 error("Unexpected failure to find reverse solution for input to output table for value %f %f %f (ICC input %f %f %f)",pp[0].v[0],pp[0].v[1],pp[0].v[2], in_v[0], in_v[1], in_v[2]);
1567 return 2;
1568 } else { /* Multiple solutions */
1569 /* Use a simple minded resolution - choose the one closest to the center */
1570 double bdist = 1e300;
1571 int bsoln = 0;
1572 DBR(("got multiple reverse solutions\n"));
1573 for (i = 0; i < nsoln; i++) {
1574 double ss;
1575
1576 DBR(("Soln %d: %s\n",i, icmPdv(p->clutTable->di, pp[i].p)))
1577 for (ss = 0.0, e = 0; e < p->clutTable->di; e++) {
1578 double tt;
1579 tt = pp[i].p[e] - p->licent[e];
1580 tt *= tt;
1581 if (tt < bdist) { /* Better solution */
1582 bdist = tt;
1583 bsoln = i;
1584 }
1585 }
1586 }
1587 #ifndef NEVER
1588 // ~~99 average them
1589 for (i = 1; i < nsoln; i++) {
1590 for (e = 0; e < p->clutTable->di; e++)
1591 pp[0].p[e] += pp[i].p[e];
1592 }
1593 for (e = 0; e < p->clutTable->di; e++)
1594 pp[0].p[e] /= (double)nsoln;
1595 bsoln = 0;
1596 #endif
1597 //printf("~1 chose %d\n",bsoln);
1598 i = bsoln;
1599 }
1600 for (e = 0; e < p->clutTable->di; e++) {
1601 /* Save solution as atractor for next one, on the basis */
1602 /* that it might have better continuity given pesudo-hilbert inversion path. */
1603 p->licent[e] = out[e] = pp[i].p[e]; /* Solution */
1604 }
1605 }
1606
1607 /* Sanitise auxiliary locus range and auxiliary value return */
1608 if (auxr != NULL || auxv != NULL) {
1609 int ee = 0;
1610 for (e = 0; e < p->clutTable->di; e++) {
1611 double v = out[e]; /* Solution */
1612 if (p->auxm[e] != 0) {
1613 if (auxr != NULL) { /* Make sure locus encloses actual value */
1614 if (auxr[2 * ee] > v)
1615 auxr[2 * ee] = v;
1616 if (auxr[2 * ee + 1] < v)
1617 auxr[2 * ee + 1] = v;
1618 }
1619 if (auxv != NULL) {
1620 auxv[ee] = v;
1621 }
1622 ee++;
1623 }
1624 }
1625 }
1626
1627 #ifdef CHECK_ILIMIT /* Do sanity checks on meeting ink limit */
1628 if (p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0) {
1629 double sum = icxLimitD(p, out);
1630 if (sum > 0.0)
1631 printf("xlut assert%s: icxLuLut_inv_clut returned outside limits by %f > tlimit %f\n",crv ? " (clip)" : "", sum, p->ink.tlimit);
1632 }
1633 #endif
1634
1635 if (clipd != NULL) {
1636 *clipd = cdist;
1637 DBR(("inv_clut_aux returning clip DE %f\n",cdist))
1638 }
1639
1640 DBR(("inv_clut_aux returning %f %f %f %f\n",out[0],out[1],out[2],out[3]))
1641 return crv;
1642 }
1643
1644 /* Do output'->input' lookup, simple version */
1645 /* Note than out[] will carry inking value if icxKrule is icxKvalue of icxKlocus */
1646 /* and that the icxKrule value will be in the input (NOT input') space. */
1647 /* Note that the ink limit will be computed after converting input' to input */
icxLuLut_inv_clut(icxLuLut * p,double * out,double * in)1648 int icxLuLut_inv_clut(icxLuLut *p, double *out, double *in) {
1649 return icxLuLut_inv_clut_aux(p, out, NULL, NULL, NULL, NULL, in);
1650 }
1651
1652 /* Given the proposed auxiliary input values in in[di], */
1653 /* and the target output' (ie. PCS') values in out[fdi], */
1654 /* return the auxiliary input (NOT input' space) values as a proportion of their */
1655 /* possible locus in locus[di]. */
1656 /* This is generally used on a source CMYK profile to convey the black intent */
1657 /* to destination CMYK profile. */
icxLuLut_clut_aux_locus(icxLuLut * p,double * locus,double * out,double * in)1658 int icxLuLut_clut_aux_locus(icxLuLut *p, double *locus, double *out, double *in) {
1659 co pp[1]; /* Room for all the solutions found */
1660 int nsoln; /* Number of solutions found */
1661 int e,f;
1662
1663 if (p->clutTable->di > p->clutTable->fdi) { /* ie. CMYK->Lab, there will be ambiguity */
1664 double min[MXDI], max[MXDI]; /* Auxiliary locus range */
1665
1666 /* Setup for auxiliary locus lookup */
1667 for (f = 0; f < p->clutTable->fdi; f++) {
1668 pp[0].v[f] = out[f]; /* Target output' (i.e. PCS) value */
1669 }
1670
1671 /* Compute auxiliary locus */
1672 nsoln = p->clutTable->rev_locus(
1673 p->clutTable, /* rspl object */
1674 p->auxm, /* Auxiliary mask */
1675 pp, /* Input target and output solutions */
1676 min, max); /* Returned locus of valid auxiliary values */
1677
1678 if (nsoln == 0) {
1679 for (e = 0; e < p->clutTable->di; e++)
1680 locus[e] = 0.0; /* Return some safe values */
1681 } else { /* Got a valid locus */
1682
1683 /* Convert the locus from input' -> input space */
1684 for (e = 0; e < p->clutTable->di; e++) {
1685 co tc;
1686 /* (Is speed more important than precision ?) */
1687 if (p->auxm[e] != 0) {
1688 tc.p[0] = min[e];
1689 p->revinputTable[e]->interp(p->revinputTable[e], &tc);
1690 min[e] = tc.v[0];
1691 tc.p[0] = max[e];
1692 p->revinputTable[e]->interp(p->revinputTable[e], &tc);
1693 max[e] = tc.v[0];
1694 }
1695 }
1696
1697 /* Figure out the proportion of the locus */
1698 for (e = 0; e < p->clutTable->di; e++) {
1699 if (p->auxm[e] != 0) {
1700 double iv = in[e];
1701 if (iv <= min[e])
1702 locus[e] = 0.0;
1703 else if (iv >= max[e])
1704 locus[e] = 1.0;
1705 else {
1706 double lpl = max[e] - min[e]; /* Locus path length */
1707 if (lpl > 1e-6)
1708 locus[e] = (iv - min[e])/lpl;
1709 else
1710 locus[e] = 0.0;
1711 }
1712 }
1713 }
1714 }
1715 } else {
1716 /* There should be no auxiliaries */
1717 for (e = 0; e < p->clutTable->di; e++)
1718 locus[e] = 0.0; /* Return some safe values */
1719 }
1720 return 0;
1721 }
1722
1723 /* Do input' -> input inverse lookup */
icxLuLut_inv_input(icxLuLut * p,double * out,double * in)1724 int icxLuLut_inv_input(icxLuLut *p, double *out, double *in) {
1725 #ifdef NEVER
1726 return ((icmLuLut *)p->plu)->inv_input((icmLuLut *)p->plu, out, in);
1727 #else
1728 int rv = 0;
1729 int i,j;
1730 int nsoln; /* Number of solutions found */
1731 co pp[MAX_INVSOLN]; /* Room for all the solutions found */
1732 // double cdir;
1733
1734 DBR(("inv_input got DEV' %f %f %f %f\n",in[0],in[1],in[2],in[3]))
1735
1736 for (i = 0; i < p->inputChan; i++) {
1737 pp[0].p[0] = p->inputClipc[i];
1738 pp[0].v[0] = in[i];
1739 // cdir = p->inputClipc[i] - in[i]; /* Clip towards output range */
1740
1741 nsoln = p->inputTable[i]->rev_interp (
1742 p->inputTable[i], /* this */
1743 RSPL_NEARCLIP, /* Clip to nearest (faster than vector) */
1744 MAX_INVSOLN, /* Maximum number of solutions allowed for */
1745 NULL, /* No auxiliary input targets */
1746 NULL, /* No LCH weight because this is 1D */
1747 // &cdir, /* Clip vector direction and length */
1748 pp); /* Input and output values */
1749
1750 if (nsoln & RSPL_DIDCLIP)
1751 rv = 1;
1752
1753 nsoln &= RSPL_NOSOLNS; /* Get number of solutions */
1754
1755 if (nsoln == 1) { /* Exactly one solution */
1756 j = 0;
1757 } else if (nsoln == 0) { /* Zero solutions. This is unexpected. */
1758 error("Unexpected failure to find reverse solution for input table");
1759 return 2;
1760 } else { /* Multiple solutions */
1761 /* Use a simple minded resolution - choose the one closest to the center */
1762 double bdist = 1e300;
1763 int bsoln = 0;
1764 /* Don't expect this - 1D luts are meant to be monotonic */
1765 warning("1D lut inversion got %d reverse solutions\n",nsoln);
1766 warning("solution 0 = %f\n",pp[0].p[0]);
1767 warning("solution 1 = %f\n",pp[1].p[0]);
1768 for (j = 0; j < nsoln; j++) {
1769 double tt;
1770 tt = pp[i].p[0] - p->inputClipc[i];
1771 tt *= tt;
1772 if (tt < bdist) { /* Better solution */
1773 bdist = tt;
1774 bsoln = j;
1775 }
1776 }
1777 j = bsoln;
1778 }
1779 out[i] = pp[j].p[0];
1780 }
1781
1782 DBR(("inv_input returning DEV %f %f %f %f\n",out[0],out[1],out[2],out[3]))
1783 return rv;
1784 #endif /* NEVER */
1785 }
1786
1787 /* Possible inverse matrix lookup */
1788 /* (Will do nothing if input is device space) */
icxLuLut_inv_matrix(icxLuLut * p,double * out,double * in)1789 int icxLuLut_inv_matrix(icxLuLut *p, double *out, double *in) {
1790 int rv = 0;
1791 rv |= ((icmLuLut *)p->plu)->inv_matrix((icmLuLut *)p->plu, out, in);
1792 return rv;
1793 }
1794
1795 /* Inverse input absolute intent conversion */
1796 /* (Will do nothing if input is device space) */
icxLuLut_inv_in_abs(icxLuLut * p,double * out,double * in)1797 int icxLuLut_inv_in_abs(icxLuLut *p, double *out, double *in) {
1798 int rv = 0;
1799 rv |= ((icmLuLut *)p->plu)->inv_in_abs((icmLuLut *)p->plu, out, in);
1800
1801 if (p->ins == icxSigJabData) {
1802 p->cam->XYZ_to_cam(p->cam, out, out);
1803 }
1804
1805 return rv;
1806 }
1807
1808 /* Overall inverse lookup */
1809 /* Note that all auxiliary values are in input (NOT input') space */
1810 static int
icxLuLut_inv_lookup(icxLuBase * pp,double * out,double * in)1811 icxLuLut_inv_lookup(
1812 icxLuBase *pp, /* This */
1813 double *out, /* Vector of output values/input auxiliary values */
1814 double *in /* Vector of input values */
1815 ) {
1816 icxLuLut *p = (icxLuLut *)pp;
1817 int rv = 0;
1818 int i;
1819 double temp[MAX_CHAN];
1820
1821 DBOL(("xiccilu: input = %s\n", icmPdv(p->outputChan, in)));
1822 if (p->mergeclut == 0) { /* Do this if it's not merger with clut */
1823 rv |= p->inv_out_abs (p, temp, in);
1824 DBOL(("xiccilu: after inv abs = %s\n", icmPdv(p->outputChan, temp)));
1825 rv |= p->inv_output (p, temp, temp);
1826 DBOL(("xiccilu: after inv out = %s\n", icmPdv(p->outputChan, temp)));
1827 } else {
1828 for (i = 0; i < p->outputChan; i++)
1829 temp[i] = in[i];
1830 }
1831 DBOL(("xiccilu: aux targ = %s\n", icmPdv(p->inputChan,out)));
1832 rv |= p->inv_clut (p, out, temp);
1833 DBOL(("xiccilu: after inv clut = %s\n", icmPdv(p->inputChan,out)));
1834 rv |= p->inv_input (p, out, out);
1835 DBOL(("xiccilu: after inv input = %s\n", icmPdv(p->inputChan,out)));
1836 rv |= p->inv_matrix (p, out, out);
1837 DBOL(("xiccilu: after inv matrix = %s\n", icmPdv(p->inputChan,out)));
1838 rv |= p->inv_in_abs (p, out, out);
1839 DBOL(("xiccilu: after inv abs = %s\n", icmPdv(p->inputChan,out)));
1840 return rv;
1841 }
1842
1843 /* - - - - - - - - - - - - - - - - - - - - - - - - - - */
1844 /* Given a nominated output PCS (ie. Absolute, Jab etc.), convert it in the bwd */
1845 /* direction into a relative XYZ or Lab PCS value */
1846 /* (This is used in generating gamut compression in B2A tables) */
icxLuLut_bwd_outpcs_relpcs(icxLuBase * pp,icColorSpaceSignature os,double * out,double * in)1847 void icxLuLut_bwd_outpcs_relpcs(
1848 icxLuBase *pp,
1849 icColorSpaceSignature os, /* Output space, XYZ or Lab */
1850 double *out, double *in) {
1851 icxLuLut *p = (icxLuLut *)pp;
1852
1853 if (p->outs == icxSigJabData) {
1854 DBS(("bwd_outpcs_relpcs: Jab in = %s\n", icmPdv(3, in)));
1855 p->cam->cam_to_XYZ(p->cam, out, in);
1856 DBS(("bwd_outpcs_relpcs: abs XYZ = %s\n", icmPdv(3, out)));
1857 /* Hack to prevent CAM02 weirdness being amplified by */
1858 /* per channel clipping. */
1859 /* Limit -Y to non-stupid values by scaling */
1860 if (out[1] < -0.1) {
1861 out[0] *= -0.1/out[1];
1862 out[2] *= -0.1/out[1];
1863 out[1] = -0.1;
1864 DBS(("bwd_outpcs_relpcs: after clipping -Y %s\n",icmPdv(p->outputChan, out)));
1865 }
1866 } else {
1867 DBS(("bwd_outpcs_relpcs: abs PCS in = %s\n", icmPdv(3, out)));
1868 icmCpy3(out, in);
1869 }
1870
1871 ((icmLuLut *)p->plu)->inv_out_abs((icmLuLut *)p->plu, out, out);
1872 DBS(("bwd_outpcs_relpcs: rel PCS = %s\n", icmPdv(3, out)));
1873
1874 if (os == icSigXYZData && p->natpcs == icSigLabData) {
1875 icmLab2XYZ(&icmD50, out, out);
1876 DBS(("bwd_outpcs_relpcs: rel XYZ = %s\n", icmPdv(3, out)));
1877 } else if (os == icSigXYZData && p->natpcs == icSigLabData) {
1878 icmXYZ2Lab(&icmD50, out, out);
1879 DBS(("bwd_outpcs_relpcs: rel Lab = %s\n", icmPdv(3, out)));
1880 }
1881 }
1882
1883 /* - - - - - - - - - - - - - - - - - - - - - - - - - - */
1884
1885 /* Return LuLut information */
icxLuLut_get_info(icxLuLut * p,icmLut ** lutp,icmXYZNumber * pcswhtp,icmXYZNumber * whitep,icmXYZNumber * blackp)1886 static void icxLuLut_get_info(
1887 icxLuLut *p, /* this */
1888 icmLut **lutp, /* Pointer to icc lut type return value */
1889 icmXYZNumber *pcswhtp, /* Pointer to profile PCS white point return value */
1890 icmXYZNumber *whitep, /* Pointer to profile absolute white point return value */
1891 icmXYZNumber *blackp /* Pointer to profile absolute black point return value */
1892 ) {
1893 ((icmLuLut *)p->plu)->get_info((icmLuLut *)p->plu, lutp, pcswhtp, whitep, blackp);
1894 }
1895
1896 /* Return the underlying Lut matrix */
1897 static void
icxLuLut_get_matrix(icxLuLut * p,double m[3][3])1898 icxLuLut_get_matrix (
1899 icxLuLut *p,
1900 double m[3][3]
1901 ) {
1902 ((icmLuLut *)p->plu)->get_matrix((icmLuLut *)p->plu, m);
1903 }
1904
1905 static void
icxLuLut_free(icxLuBase * pp)1906 icxLuLut_free(
1907 icxLuBase *pp
1908 ) {
1909 icxLuLut *p = (icxLuLut *)pp;
1910 int i;
1911
1912 for (i = 0; i < p->inputChan; i++) {
1913 if (p->inputTable[i] != NULL)
1914 p->inputTable[i]->del(p->inputTable[i]);
1915 if (p->revinputTable[i] != NULL)
1916 p->revinputTable[i]->del(p->revinputTable[i]);
1917 }
1918
1919 if (p->clutTable != NULL)
1920 p->clutTable->del(p->clutTable);
1921
1922 if (p->cclutTable != NULL)
1923 p->cclutTable->del(p->cclutTable);
1924
1925 for (i = 0; i < p->outputChan; i++) {
1926 if (p->outputTable[i] != NULL)
1927 p->outputTable[i]->del(p->outputTable[i]);
1928 }
1929
1930 if (p->plu != NULL)
1931 p->plu->del(p->plu);
1932
1933 if (p->cam != NULL)
1934 p->cam->del(p->cam);
1935
1936 if (p->absxyzlu != NULL)
1937 p->absxyzlu->del(p->absxyzlu);
1938
1939 free(p);
1940 }
1941
1942 /* - - - - - - - - - - - - - - - - - - - - - - - - - - */
1943
1944 static gamut *icxLuLutGamut(icxLuBase *plu, double detail);
1945 static icxCuspMap *icxLuLutCuspMap(icxLuBase *plu, int res);
1946
1947 /* Do the basic icxLuLut creation and initialisation */
1948 static icxLuLut *
alloc_icxLuLut(xicc * xicp,icmLuBase * plu,int flags)1949 alloc_icxLuLut(
1950 xicc *xicp,
1951 icmLuBase *plu, /* Pointer to Lu we are expanding (ours) */
1952 int flags /* clip, merge flags */
1953 ) {
1954 icxLuLut *p; /* Object being created */
1955 icmLuLut *luluto = (icmLuLut *)plu; /* Lookup Lut type object */
1956
1957 if ((p = (icxLuLut *) calloc(1,sizeof(icxLuLut))) == NULL)
1958 return NULL;
1959
1960 p->pp = xicp;
1961 p->plu = plu;
1962 p->del = icxLuLut_free;
1963 p->lutspaces = icxLutSpaces;
1964 p->spaces = icxLuSpaces;
1965 p->get_native_ranges = icxLu_get_native_ranges;
1966 p->get_ranges = icxLu_get_ranges;
1967 p->efv_wh_bk_points = icxLuEfv_wh_bk_points;
1968 p->get_gamut = icxLuLutGamut;
1969 p->get_cuspmap = icxLuLutCuspMap;
1970 p->fwd_relpcs_outpcs = icxLuLut_fwd_relpcs_outpcs;
1971 p->bwd_outpcs_relpcs = icxLuLut_bwd_outpcs_relpcs;
1972 p->nearclip = 0; /* Set flag defaults */
1973 p->mergeclut = 0;
1974 p->noisluts = 0;
1975 p->noipluts = 0;
1976 p->nooluts = 0;
1977 p->intsep = 0;
1978
1979 p->lookup = icxLuLut_lookup;
1980 p->in_abs = icxLuLut_in_abs;
1981 p->matrix = icxLuLut_matrix;
1982 p->input = icxLuLut_input;
1983 p->clut = icxLuLut_clut;
1984 p->clut_aux = icxLuLut_clut_aux;
1985 p->output = icxLuLut_output;
1986 p->out_abs = icxLuLut_out_abs;
1987
1988 p->inv_lookup = icxLuLut_inv_lookup;
1989 p->inv_in_abs = icxLuLut_inv_in_abs;
1990 p->inv_matrix = icxLuLut_inv_matrix;
1991 p->inv_input = icxLuLut_inv_input;
1992 p->inv_clut = icxLuLut_inv_clut;
1993 p->inv_clut_aux = icxLuLut_inv_clut_aux;
1994 p->inv_output = icxLuLut_inv_output;
1995 p->inv_out_abs = icxLuLut_inv_out_abs;
1996
1997 p->clut_locus = icxLuLut_clut_aux_locus;
1998
1999 p->get_info = icxLuLut_get_info;
2000 p->get_matrix = icxLuLut_get_matrix;
2001
2002 /* Setup all the rspl analogs of the icc Lut */
2003 /* NOTE: We assume that none of this relies on the flag settings, */
2004 /* since they will be set on our return. */
2005
2006 /* Get details of underlying, native icc color space */
2007 p->plu->lutspaces(p->plu, &p->natis, NULL, &p->natos, NULL, &p->natpcs);
2008
2009 /* Get other details of conversion */
2010 p->plu->spaces(p->plu, NULL, &p->inputChan, NULL, &p->outputChan, NULL, NULL, NULL, NULL, NULL);
2011
2012 /* Remember flags */
2013 p->flags = flags;
2014
2015 /* Sanity check */
2016 if (p->inputChan > MXDI) {
2017 sprintf(p->pp->err,"xicc can only handle input channels of %d or less",MXDI);
2018 p->inputChan = MXDI; /* Avoid going outside array bounds */
2019 p->pp->errc = 1;
2020 p->del((icxLuBase *)p);
2021 return NULL;
2022 }
2023 if (p->outputChan > MXDO) {
2024 sprintf(p->pp->err,"xicc can only handle output channels of %d or less",MXDO);
2025 p->outputChan = MXDO; /* Avoid going outside array bounds */
2026 p->pp->errc = 1;
2027 p->del((icxLuBase *)p);
2028 return NULL;
2029 }
2030
2031 /* Get pointer to icmLut */
2032 luluto->get_info(luluto, &p->lut, NULL, NULL, NULL);
2033
2034 return p;
2035 }
2036
2037 /* Initialise the clut ink limiting and black */
2038 /* generation information. */
2039 /* return 0 or error code */
2040 static int
setup_ink_icxLuLut(icxLuLut * p,icxInk * ink,int setLminmax)2041 setup_ink_icxLuLut(
2042 icxLuLut *p, /* Object being initialised */
2043 icxInk *ink, /* inking details (NULL for default) */
2044 int setLminmax /* Figure the L locus for inking rule */
2045 ) {
2046 int devchan = p->func == icmFwd ? p->inputChan : p->outputChan;
2047
2048 if (ink) {
2049 p->ink = *ink; /* Copy the structure */
2050 } else {
2051 p->ink.tlimit = 3.0; /* default ink limit of 300% */
2052 p->ink.klimit = -1.0; /* default no black limit */
2053 p->ink.KonlyLmin = 0; /* default don't use K only black as Locus Lmin */
2054 p->ink.k_rule = icxKluma5; /* default K generation rule */
2055 p->ink.c.Ksmth = ICXINKDEFSMTH; /* Default smoothing */
2056 p->ink.c.Kskew = ICXINKDEFSKEW; /* Default shape skew */
2057 p->ink.c.Kstle = 0.0; /* Min K at white end */
2058 p->ink.c.Kstpo = 0.0; /* Start of transition is at white */
2059 p->ink.c.Kenle = 1.0; /* Max K at black end */
2060 p->ink.c.Kenpo = 1.0; /* End transition at black */
2061 p->ink.c.Kshap = 1.0; /* Linear transition */
2062 }
2063
2064 /* Normalise total and black ink limits */
2065 if (p->ink.tlimit <= 1e-4 || p->ink.tlimit >= (double)devchan)
2066 p->ink.tlimit = -1.0; /* Turn ink limit off if not effective */
2067 if (devchan < 4 || p->ink.klimit < 0.0 || p->ink.klimit >= 1.0)
2068 p->ink.klimit = -1.0; /* Turn black limit off if not effective */
2069
2070 /* Set the ink limit information for any reverse interpolation. */
2071 /* Calling this will clear the reverse interpolaton cache. */
2072 p->clutTable->rev_set_limit(
2073 p->clutTable, /* this */
2074 p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0 ? icxLimitD_void : NULL,
2075 /* Optional input space limit() function. */
2076 /* Function should evaluate in[0..di-1], and return */
2077 /* number that is not to exceed 0.0. NULL if not used. */
2078 (void *)p, /* Context passed to limit() */
2079 0.0 /* Value that limit() is not to exceed */
2080 );
2081
2082 /* Duplicate in the CAM clip rspl if it exists */
2083 if (p->cclutTable != NULL) {
2084 p->cclutTable->rev_set_limit(
2085 p->cclutTable, /* this */
2086 p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0 ? icxLimitD_void : NULL,
2087 /* Optional input space limit() function. */
2088 /* Function should evaluate in[0..di-1], and return */
2089 /* number that is not to exceed 0.0. NULL if not used. */
2090 (void *)p, /* Context passed to limit() */
2091 0.0 /* Value that limit() is not to exceed */
2092 );
2093 }
2094
2095 /* Figure Lmin and Lmax for icxKluma5 curve basis */
2096 if (setLminmax
2097 && p->clutTable->di > p->clutTable->fdi) { /* If K generation makes sense */
2098 double wh[3], bk[3], kk[3];
2099 int mergeclut; /* Save/restore mergeclut value */
2100
2101 /* Get white/black in effective xlu PCS space */
2102 p->efv_wh_bk_points((icxLuBase *)p, wh, bk, kk);
2103
2104 /* Convert from effective PCS (ie. Jab) to native XYZ or Lab PCS */
2105 mergeclut = p->mergeclut; /* Hack to be able to use inv_out_abs() */
2106 p->mergeclut = 0; /* if mergeclut is active. */
2107 icxLuLut_inv_out_abs(p, wh, wh);
2108 icxLuLut_inv_out_abs(p, bk, bk);
2109 icxLuLut_inv_out_abs(p, kk, kk);
2110 p->mergeclut = mergeclut; /* Restore */
2111
2112 /* Convert to Lab PCS */
2113 if (p->natos == icSigXYZData) { /* Always do K rule in L space */
2114 icmXYZ2Lab(&icmD50, wh, wh);
2115 icmXYZ2Lab(&icmD50, bk, bk);
2116 icmXYZ2Lab(&icmD50, kk, kk);
2117 }
2118 p->Lmax = 0.01 * wh[0];
2119 if (p->ink.KonlyLmin != 0)
2120 p->Lmin = 0.01 * kk[0];
2121 else
2122 p->Lmin = 0.01 * bk[0];
2123 } else {
2124
2125 /* Some sane defaults */
2126 p->Lmax = 1.0;
2127 p->Lmin = 0.0;
2128 }
2129
2130 return 0;
2131 }
2132
2133
2134 /* Initialise the clut clipping information, ink limiting */
2135 /* and auxiliary parameter settings for all the rspl. */
2136 /* return 0 or error code */
2137 static int
setup_clip_icxLuLut(icxLuLut * p)2138 setup_clip_icxLuLut(
2139 icxLuLut *p /* Object being initialised */
2140 ) {
2141 double tmin[MXDIDO], tmax[MXDIDO];
2142 int i;
2143
2144 /* Setup for inversion of multi-dim clut */
2145
2146 p->kch = -1; /* Not known yet */
2147
2148 /* Set auxiliaries */
2149 for (i = 0; i < p->inputChan; i++)
2150 p->auxm[i] = 0;
2151
2152 if (p->inputChan > p->outputChan) {
2153 switch(p->natis) {
2154 case icSigCmykData:
2155 /* Should fix icm/xicc to remember K channel */
2156 p->auxm[3] = 1; /* K is the auxiliary channel */
2157 break;
2158 default:
2159 if (p->kch >= 0) /* It's been discovered */
2160 p->auxm[p->kch] = 1;
2161 else {
2162 p->pp->errc = 2;
2163 sprintf(p->pp->err,"Unknown colorspace %s when setting auxliaries",
2164 icm2str(icmColorSpaceSignature, p->natis));
2165 return p->pp->errc;
2166 }
2167 break;
2168 }
2169 }
2170
2171 // p->auxlinf = NULL; /* Input space auxiliary linearisation function - not implemented */
2172 // p->auxlinf_ctx = NULL; /* Opaque context for auxliliary linearisation function */
2173
2174 /* Aproximate center of clut input gamut - used for */
2175 /* resolving multiple reverse solutions. */
2176 p->clutTable->get_in_range(p->clutTable, tmin, tmax);
2177 for (i = 0; i < p->clutTable->di; i++) {
2178 p->licent[i] = p->icent[i] = (tmin[i] + tmax[i])/2.0;
2179 }
2180
2181 /* Compute clip setup information relating to clut output gamut. */
2182 if (p->nearclip != 0 /* Near clip requested */
2183 || p->inputChan == 1) { /* or vector clip won't work */
2184 p->clip.nearclip = 1;
2185
2186 } else { /* Vector clip */
2187 icColorSpaceSignature clutos = p->natos;
2188
2189 p->clip.nearclip = 0;
2190 p->clip.LabLike = 0;
2191 p->clip.fdi = p->clutTable->fdi;
2192
2193 switch(clutos) {
2194 case icxSigJabData:
2195 case icSigLabData: {
2196
2197 p->clip.LabLike = 1;
2198
2199 /* Create a CuspMap to point vectors towards */
2200 /* (Don't make it too fine, or there will be dips) */
2201 p->clip.cm = p->get_cuspmap((icxLuBase *)p, 30);
2202 break;
2203 }
2204 case icSigXYZData:
2205 // ~~~~~~1 need to add this.
2206 warning("xlut.c: setup_clip_icxLuLut() icSigXYZData case not implemented!");
2207 /* Fall through */
2208
2209 default:
2210 /* Do a crude approximation, that may not work. */
2211 p->clutTable->get_out_range(p->clutTable, tmin, tmax);
2212 for (i = 0; i < p->clutTable->fdi; i++)
2213 p->clip.ocent[i] = (tmin[i] + tmax[i])/2.0;
2214
2215 break;
2216 }
2217 }
2218 return 0;
2219 }
2220
2221 /* Function to pass to rspl to set secondary input/output transfer functions */
2222 static void
icxLuLut_inout_func(void * pp,double * out,double * in)2223 icxLuLut_inout_func(
2224 void *pp, /* icxLuLut */
2225 double *out, /* output value */
2226 double *in /* inut value */
2227 ) {
2228 icxLuLut *p = (icxLuLut *)pp; /* this */
2229 icmLuLut *luluto = (icmLuLut *)p->plu; /* Get icmLuLut object */
2230 double tin[MAX_CHAN];
2231 double tout[MAX_CHAN];
2232 int i;
2233
2234 if (p->iol_out == 0) { /* fwd input */
2235 #ifdef INK_LIMIT_TEST
2236 tout[p->iol_ch] = in[0];
2237 #else
2238 for (i = 0; i < p->inputChan; i++)
2239 tin[i] = 0.0;
2240 tin[p->iol_ch] = in[0];
2241 luluto->input(luluto, tout, tin);
2242 #endif
2243 } else if (p->iol_out == 1) { /* fwd output */
2244 for (i = 0; i < p->outputChan; i++)
2245 tin[i] = 0.0;
2246 tin[p->iol_ch] = in[0];
2247 luluto->output(luluto, tout, tin);
2248 } else { /* bwd input */
2249 #ifdef INK_LIMIT_TEST
2250 tout[p->iol_ch] = in[0];
2251 #else
2252 for (i = 0; i < p->inputChan; i++)
2253 tin[i] = 0.0;
2254 tin[p->iol_ch] = in[0];
2255 luluto->inv_input(luluto, tout, tin);
2256 /* This won't be valid if matrix is used or there is a PCS conversion !!! */
2257 /* Shouldn't be a problem because this is only used for inverting dev->PCS ? */
2258 luluto->inv_in_abs(luluto, tout, tout);
2259 #endif
2260 }
2261 out[0] = tout[p->iol_ch];
2262 }
2263
2264 /* Function to pass to rspl to set clut up, when mergeclut is set */
2265 static void
icxLuLut_clut_merge_func(void * pp,double * out,double * in)2266 icxLuLut_clut_merge_func(
2267 void *pp, /* icxLuLut */
2268 double *out, /* output value */
2269 double *in /* input value */
2270 ) {
2271 icxLuLut *p = (icxLuLut *)pp; /* this */
2272 icmLuLut *luluto = (icmLuLut *)p->plu; /* Get icmLuLut object */
2273
2274 luluto->clut(luluto, out, in);
2275 luluto->output(luluto, out, out);
2276 luluto->out_abs(luluto, out, out);
2277
2278 if (p->outs == icxSigJabData) {
2279 p->cam->XYZ_to_cam(p->cam, out, out);
2280 }
2281 }
2282
2283 /* Implimenation of Lut create from icc. */
2284 /* Note that xicc_get_luobj() will have set the pcsor & */
2285 /* intent to consistent values if Jab and/or icxAppearance */
2286 /* has been requested. */
2287 /* It will also have created the underlying icm lookup object */
2288 /* that is used to create and implement the icx one. The icm */
2289 /* will be used to translate from native to effective PCS, unless */
2290 /* the effective PCS is Jab, in which case the icm will be set to */
2291 /* have an effective PCS of XYZ. Since native<->effecive PCS conversion */
2292 /* is done at the to/from_abs() stage, none of it affects the individual */
2293 /* conversion steps, which will all talk the native PCS (unless merged). */
2294 static icxLuBase *
new_icxLuLut(xicc * xicp,int flags,icmLuBase * plu,icmLookupFunc func,icRenderingIntent intent,icColorSpaceSignature pcsor,icxViewCond * vc,icxInk * ink)2295 new_icxLuLut(
2296 xicc *xicp,
2297 int flags, /* clip, merge flags */
2298 icmLuBase *plu, /* Pointer to Lu we are expanding (ours) */
2299 icmLookupFunc func, /* Functionality requested */
2300 icRenderingIntent intent, /* Rendering intent */
2301 icColorSpaceSignature pcsor, /* PCS override (0 = def) */
2302 icxViewCond *vc, /* Viewing Condition (NULL if not using CAM) */
2303 icxInk *ink /* inking details (NULL for default) */
2304 ) {
2305 icxLuLut *p; /* Object being created */
2306 icmLuLut *luluto = (icmLuLut *)plu; /* Lookup Lut type object */
2307 icmLookupFunc fnc;
2308
2309 int i;
2310
2311 /* Do basic creation and initialisation */
2312 if ((p = alloc_icxLuLut(xicp, plu, flags)) == NULL)
2313 return NULL;
2314
2315 p->func = func;
2316
2317 /* Set LuLut "use" specific creation flags: */
2318 if (flags & ICX_CLIP_NEAREST)
2319 p->nearclip = 1;
2320
2321 if (flags & ICX_MERGE_CLUT)
2322 p->mergeclut = 1;
2323
2324 if (flags & ICX_FAST_SETUP)
2325 p->fastsetup = 1;
2326
2327 /* We're only implementing this under specific conditions. */
2328 if (flags & ICX_CAM_CLIP
2329 && func == icmFwd
2330 && !(p->mergeclut != 0 && pcsor == icxSigJabData)) /* Don't need camclip if merged Jab */
2331 p->camclip = 1;
2332
2333 if (flags & ICX_INT_SEPARATE) {
2334 fprintf(stderr,"~1 Internal optimised 4D separations not yet implemented!\n");
2335 p->intsep = 1;
2336 }
2337
2338 /* Init the CAM model if it will be used */
2339 if (pcsor == icxSigJabData || p->camclip) {
2340 if (vc != NULL) /* One has been provided */
2341 p->vc = *vc; /* Copy the structure */
2342 else
2343 xicc_enum_viewcond(xicp, &p->vc, -1, NULL, 0, NULL); /* Use a default */
2344 p->cam = new_icxcam(cam_default);
2345 p->cam->set_view(p->cam, p->vc.Ev, p->vc.Wxyz, p->vc.La, p->vc.Yb, p->vc.Lv,
2346 p->vc.Yf, p->vc.Yg, p->vc.Gxyz, XICC_USE_HK, p->vc.hkscale);
2347 } else {
2348 p->cam = NULL;
2349 }
2350
2351 /* Remember the effective intent */
2352 p->intent = intent;
2353
2354 /* Get the effective spaces of underlying icm */
2355 plu->spaces(plu, &p->ins, NULL, &p->outs, NULL, NULL, NULL, &fnc, &p->pcs, NULL);
2356
2357 /* Override with pcsor */
2358 /* We assume that any profile that has a CIE color as a "device" color */
2359 /* intends it to stay that way, and not be overridden. */
2360 if (pcsor == icxSigJabData) {
2361 p->pcs = pcsor;
2362
2363 if (xicp->pp->header->deviceClass == icSigAbstractClass) {
2364 p->ins = pcsor;
2365 p->outs = pcsor;
2366
2367 } else if (xicp->pp->header->deviceClass != icSigLinkClass) {
2368 if (func == icmBwd || func == icmGamut || func == icmPreview)
2369 p->ins = pcsor;
2370 if (func == icmFwd || func == icmPreview)
2371 p->outs = pcsor;
2372 }
2373 }
2374
2375 /* In general the native and effective ranges of the icx will be the same as the */
2376 /* underlying icm lookup object. */
2377 p->plu->get_lutranges(p->plu, p->ninmin, p->ninmax, p->noutmin, p->noutmax);
2378 p->plu->get_ranges(p->plu, p->inmin, p->inmax, p->outmin, p->outmax);
2379
2380 /* If we have a Jab PCS override, reflect this in the effective icx range. */
2381 /* Note that the ab ranges are nominal. They will exceed this range */
2382 /* for colors representable in L*a*b* PCS */
2383 if (p->ins == icxSigJabData) {
2384 p->inmin[0] = 0.0; p->inmax[0] = 100.0;
2385 p->inmin[1] = -128.0; p->inmax[1] = 128.0;
2386 p->inmin[2] = -128.0; p->inmax[2] = 128.0;
2387 } else if (p->outs == icxSigJabData) {
2388 p->outmin[0] = 0.0; p->outmax[0] = 100.0;
2389 p->outmin[1] = -128.0; p->outmax[1] = 128.0;
2390 p->outmin[2] = -128.0; p->outmax[2] = 128.0;
2391 }
2392
2393 /* If we have a merged clut, reflect this in the icx native PCS range. */
2394 /* Merging merges output processing (irrespective of whether we are using */
2395 /* the forward or backward cluts) */
2396 if (p->mergeclut != 0) {
2397 int i;
2398 for (i = 0; i < p->outputChan; i++) {
2399 p->noutmin[i] = p->outmin[i];
2400 p->noutmax[i] = p->outmax[i];
2401 }
2402 }
2403
2404 /* ------------------------------- */
2405 /* Create rspl based input lookups */
2406 for (i = 0; i < p->inputChan; i++) {
2407 if ((p->inputTable[i] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) {
2408 p->pp->errc = 2;
2409 sprintf(p->pp->err,"Creation of input table rspl failed");
2410 p->del((icxLuBase *)p);
2411 return NULL;
2412 }
2413 p->iol_out = 0; /* Input lookup */
2414 p->iol_ch = i; /* Chanel */
2415 p->inputTable[i]->set_rspl(p->inputTable[i], RSPL_NOFLAGS,
2416 (void *)p, icxLuLut_inout_func,
2417 &p->ninmin[i], &p->ninmax[i], (int *)&p->lut->inputEnt, &p->ninmin[i], &p->ninmax[i]);
2418 }
2419
2420 /* Setup center clip target for inversion */
2421 for (i = 0; i < p->inputChan; i++) {
2422 p->inputClipc[i] = (p->ninmin[i] + p->ninmax[i])/2.0;
2423 }
2424
2425 /* Create rspl based reverse input lookups used in ink limit and locus range functions. */
2426 for (i = 0; i < p->inputChan; i++) {
2427 int gres = p->inputTable[i]->g.mres; /* Same res as curve being inverted */
2428 if (gres < 256)
2429 gres = 256;
2430 if ((p->revinputTable[i] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) {
2431 p->pp->errc = 2;
2432 sprintf(p->pp->err,"Creation of reverse input table rspl failed");
2433 p->del((icxLuBase *)p);
2434 return NULL;
2435 }
2436 p->iol_out = 2; /* Input lookup */
2437 p->iol_ch = i; /* Chanel */
2438 p->revinputTable[i]->set_rspl(p->revinputTable[i], RSPL_NOFLAGS,
2439 (void *)p, icxLuLut_inout_func,
2440 &p->ninmin[i], &p->ninmax[i], &gres, &p->ninmin[i], &p->ninmax[i]);
2441 }
2442
2443 /* ------------------------------- */
2444 {
2445 int gres[MXDI];
2446 int xflags = 0;
2447
2448 for (i = 0; i < p->inputChan; i++)
2449 gres[i] = p->lut->clutPoints;
2450
2451 #ifdef FASTREVSETUP_NON_CAM
2452 /* Don't fill in nnrev array if we aren't going to use it */
2453 if (p->camclip && p->nearclip)
2454 xflags = RSPL_FASTREVSETUP;
2455 #endif
2456
2457 /* Create rspl based multi-dim table */
2458 if ((p->clutTable = new_rspl((p->fastsetup ? RSPL_FASTREVSETUP : RSPL_NOFLAGS)
2459 | (flags & ICX_VERBOSE ? RSPL_VERBOSE : RSPL_NOFLAGS)
2460 | xflags,
2461 p->inputChan, p->outputChan)) == NULL) {
2462 p->pp->errc = 2;
2463 sprintf(p->pp->err,"Creation of clut table rspl failed");
2464 p->del((icxLuBase *)p);
2465 return NULL;
2466 }
2467
2468 if (p->mergeclut == 0) { /* Do this if it's not merged with clut, */
2469 p->clutTable->set_rspl(p->clutTable, RSPL_NOFLAGS,
2470 (void *)luluto, (void (*)(void *, double *, double *))luluto->clut,
2471 p->ninmin, p->ninmax, gres, p->noutmin, p->noutmax);
2472
2473 } else { /* If mergeclut */
2474 p->clutTable->set_rspl(p->clutTable, RSPL_NOFLAGS,
2475 (void *)p, icxLuLut_clut_merge_func,
2476 p->ninmin, p->ninmax, gres, p->noutmin, p->noutmax);
2477
2478 }
2479
2480 #ifdef USELCHWEIGHT
2481 /* If we are not doing camclip, but our output is an Lab like space, */
2482 /* then apply lchw weighting anyway. */
2483 if (!p->camclip && (p->outs == icSigLabData || p->outs == icxSigJabData)) {
2484 double lchw[MXRO] = { JCCWEIGHT, CCCWEIGHT, HCCWEIGHT };
2485
2486 /* Set the Nearest Neighbor clipping Weighting */
2487 p->clutTable->rev_set_lchw(p->clutTable, lchw);
2488 }
2489 #endif /* USELCHWEIGHT */
2490
2491 /* clut clipping is setup separately */
2492 }
2493
2494 /* ------------------------------- */
2495 /* Create rspl based output lookups */
2496 for (i = 0; i < p->outputChan; i++) {
2497 if ((p->outputTable[i] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) {
2498 p->pp->errc = 2;
2499 sprintf(p->pp->err,"Creation of output table rspl failed");
2500 p->del((icxLuBase *)p);
2501 return NULL;
2502 }
2503 p->iol_out = 1; /* Output lookup */
2504 p->iol_ch = i; /* Chanel */
2505 p->outputTable[i]->set_rspl(p->outputTable[i], RSPL_NOFLAGS,
2506 (void *)p, icxLuLut_inout_func,
2507 &p->noutmin[i], &p->noutmax[i], (int *)&p->lut->outputEnt, &p->noutmin[i], &p->noutmax[i]);
2508 }
2509
2510 /* Setup center clip target for inversion */
2511 for (i = 0; i < p->outputChan; i++) {
2512 p->outputClipc[i] = (p->noutmin[i] + p->noutmax[i])/2.0;
2513 }
2514
2515 /* ------------------------------- */
2516
2517 /* Setup all the clipping, ink limiting and auxiliary stuff, */
2518 /* in case a reverse call is used. Only do this if we know */
2519 /* the reverse stuff isn't going to fail due to channel limits. */
2520 if (fnc != icmGamut && fnc != icmPreview
2521 && p->clutTable->within_restrictedsize(p->clutTable)) {
2522
2523 if (setup_ink_icxLuLut(p, ink, 1) != 0) {
2524 p->del((icxLuBase *)p);
2525 return NULL;
2526 }
2527
2528 if (setup_clip_icxLuLut(p) != 0) {
2529 p->del((icxLuBase *)p);
2530 return NULL;
2531 }
2532 }
2533
2534 return (icxLuBase *)p;
2535 }
2536
2537
2538 /* Function to pass to rspl to set clut up, when camclip is going to be used. */
2539 /* We use the temporary icm fwd absolute xyz lookup, then convert to CAM Jab. */
2540 static void
icxLuLut_clut_camclip_func(void * pp,double * out,double * in)2541 icxLuLut_clut_camclip_func(
2542 void *pp, /* icxLuLut */
2543 double *out, /* output value */
2544 double *in /* inut value */
2545 ) {
2546 icxLuLut *p = (icxLuLut *)pp; /* this */
2547 icmLuLut *luluto = (icmLuLut *)p->absxyzlu;
2548
2549 luluto->clut(luluto, out, in);
2550 luluto->output(luluto, out, out);
2551 luluto->out_abs(luluto, out, out);
2552 p->cam->XYZ_to_cam(p->cam, out, out);
2553 }
2554
2555 /* Initialise the additional CAM space clut rspl, used to compute */
2556 /* reverse lookup CAM clipping results when the camclip flag is set. */
2557 /* We weight the CAM nn clipping, to give a more L* and H* preserving clip direction. */
2558 /* Return error code. */
2559 /* (We are assuming nearest clipping - we aren't setting up properly for */
2560 /* vector clipping) */
2561 static int
icxLuLut_init_clut_camclip(icxLuLut * p)2562 icxLuLut_init_clut_camclip(
2563 icxLuLut *p) {
2564 int e, gres[MXDI];
2565 double lchw[MXRO] = { JCCWEIGHT, CCCWEIGHT, HCCWEIGHT };
2566
2567 /* Setup so clut contains transform to CAM Jab */
2568 /* (camclip is only used in fwd or invfwd direction lookup) */
2569 double cmin[3], cmax[3];
2570 cmin[0] = 0.0; cmax[0] = 100.0; /* Nominal Jab output ranges */
2571 cmin[1] = -128.0; cmax[1] = 128.0;
2572 cmin[2] = -128.0; cmax[2] = 128.0;
2573
2574 /* Get icm lookup we need for seting up and using CAM icx clut */
2575 if ((p->absxyzlu = p->pp->pp->get_luobj(p->pp->pp, icmFwd, icAbsoluteColorimetric,
2576 icSigXYZData, icmLuOrdNorm)) == NULL) {
2577 p->pp->errc = p->pp->pp->errc; /* Copy error to xicc */
2578 strcpy(p->pp->err, p->pp->pp->err);
2579 return p->pp->errc;
2580 }
2581
2582 /* Create CAM rspl based multi-dim table */
2583 if ((p->cclutTable = new_rspl((p->fastsetup ? RSPL_FASTREVSETUP : RSPL_NOFLAGS)
2584 | (p->flags & ICX_VERBOSE ? RSPL_VERBOSE : RSPL_NOFLAGS),
2585 p->inputChan, p->outputChan)) == NULL) {
2586 p->pp->errc = 2;
2587 sprintf(p->pp->err,"Creation of clut table rspl failed");
2588 return p->pp->errc;
2589 }
2590
2591 #ifdef USELCHWEIGHT
2592 /* Set the Nearest Neighbor clipping Weighting */
2593 p->cclutTable->rev_set_lchw(p->cclutTable, lchw);
2594 #endif /* USELCHWEIGHT */
2595
2596 for (e = 0; e < p->inputChan; e++)
2597 gres[e] = p->lut->clutPoints;
2598
2599 /* Setup our special CAM space rspl */
2600 p->cclutTable->set_rspl(p->cclutTable, RSPL_NOFLAGS, (void *)p,
2601 icxLuLut_clut_camclip_func,
2602 p->ninmin, p->ninmax, gres, cmin, cmax);
2603
2604 /* Duplicate the ink limit information for any reverse interpolation. */
2605 p->cclutTable->rev_set_limit(
2606 p->cclutTable, /* this */
2607 p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0 ? icxLimitD_void : NULL,
2608 /* Optional input space limit() function. */
2609 /* Function should evaluate in[0..di-1], and return */
2610 /* number that is not to exceed 0.0. NULL if not used. */
2611 (void *)p, /* Context passed to limit() */
2612 0.0 /* Value that limit() is not to exceed */
2613 );
2614 return 0;
2615 }
2616
2617 /* ========================================================== */
2618 /* xicc creation code */
2619 /* ========================================================== */
2620
2621 /* Callback for computing delta E squared for two output (PCS) */
2622 /* values, passed as a callback to xfit */
xfit_to_de2(void * cntx,double * in1,double * in2)2623 static double xfit_to_de2(void *cntx, double *in1, double *in2) {
2624 icxLuLut *p = (icxLuLut *)cntx;
2625 double rv;
2626
2627 if (p->pcs == icSigLabData) {
2628 #ifdef USE_CIE94_DE
2629 rv = icmCIE94sq(in1, in2);
2630 #else
2631 rv = icmLabDEsq(in1, in2);
2632 #endif
2633 } else {
2634 double lab1[3], lab2[3];
2635 icmXYZ2Lab(&icmD50, lab1, in1);
2636 icmXYZ2Lab(&icmD50, lab2, in2);
2637 #ifdef USE_CIE94_DE
2638 rv = icmCIE94sq(lab1, lab2);
2639 #else
2640 rv = icmLabDEsq(lab1, lab2);
2641 #endif
2642 }
2643 return rv;
2644 }
2645
2646 /* Same as above plus partial derivatives */
xfit_to_dde2(void * cntx,double dout[2][MXDIDO],double * in1,double * in2)2647 static double xfit_to_dde2(void *cntx, double dout[2][MXDIDO], double *in1, double *in2) {
2648 icxLuLut *p = (icxLuLut *)cntx;
2649 double rv;
2650
2651 if (p->pcs == icSigLabData) {
2652 int j,k;
2653 double tdout[2][3];
2654 #ifdef USE_CIE94_DE
2655 rv = icxdCIE94sq(tdout, in1, in2);
2656 #else
2657 rv = icxdLabDEsq(tdout, in1, in2);
2658 #endif
2659 for (k = 0; k < 2; k++) {
2660 for (j = 0; j < 3; j++)
2661 dout[k][j] = tdout[k][j];
2662 }
2663 } else {
2664 double lab1[3], lab2[3];
2665 double dout12[2][3][3];
2666 double tdout[2][3];
2667 int i,j,k;
2668
2669 icxdXYZ2Lab(&icmD50, lab1, dout12[0], in1);
2670 icxdXYZ2Lab(&icmD50, lab2, dout12[1], in2);
2671 #ifdef USE_CIE94_DE
2672 rv = icxdCIE94sq(tdout, lab1, lab2);
2673 #else
2674 rv = icxdLabDEsq(tdout, lab1, lab2);
2675 #endif
2676 /* Compute partial derivative (is this correct ??) */
2677 for (k = 0; k < 2; k++) {
2678 for (j = 0; j < 3; j++) {
2679 dout[k][j] = 0.0;
2680 for (i = 0; i < 3; i++) {
2681 dout[k][j] += tdout[k][i] * dout12[k][i][j];
2682 }
2683 }
2684 }
2685 }
2686 return rv;
2687 }
2688
2689 #ifdef NEVER
2690 /* Check partial derivative function within xfit_to_dde2() */
2691
_xfit_to_dde2(void * cntx,double dout[2][MXDIDO],double * in1,double * in2)2692 static double _xfit_to_dde2(void *cntx, double dout[2][MXDIDO], double *in1, double *in2) {
2693 icxLuLut *pp = (icxLuLut *)cntx;
2694 int k, i;
2695 double rv, drv;
2696 double trv;
2697
2698 rv = xfit_to_de2(cntx, in1, in2);
2699 drv = xfit_to_dde2(cntx, dout, in1, in2);
2700
2701 if (fabs(rv - drv) > 1e-6)
2702 printf("######## DDE2: RV MISMATCH is %f should be %f ########\n",rv,drv);
2703
2704 /* Check each parameter delta */
2705 for (k = 0; k < 2; k++) {
2706 for (i = 0; i < 3; i++) {
2707 double *in;
2708 double del;
2709
2710 if (k == 0)
2711 in = in1;
2712 else
2713 in = in2;
2714
2715 in[i] += 1e-9;
2716 trv = xfit_to_de2(cntx, in1, in2);
2717 in[i] -= 1e-9;
2718
2719 /* Check that del is correct */
2720 del = (trv - rv)/1e-9;
2721 if (fabs(dout[k][i] - del) > 0.04) {
2722 printf("######## DDE2: EXCESSIVE at in[%d][%d] is %f should be %f ########\n",k,i,dout[k][i],del);
2723 }
2724 }
2725 }
2726 return rv;
2727 }
2728
2729 #define xfit_to_dde2 _xfit_to_dde2
2730
2731 #endif
2732
2733 /* Context for rspl setting input and output curves */
2734 typedef struct {
2735 int iix;
2736 int oix;
2737 xfit *xf; /* Optimisation structure */
2738 } curvectx;
2739
2740 /* Function to pass to rspl to set input and output */
2741 /* transfer function for xicc lookup function */
2742 static void
set_linfunc(void * cc,double * out,double * in)2743 set_linfunc(
2744 void *cc, /* curvectx structure */
2745 double *out, /* Device output value */
2746 double *in /* Device input value */
2747 ) {
2748 curvectx *c = (curvectx *)cc; /* this */
2749 xfit *p = c->xf;
2750
2751 if (c->iix >= 0) { /* Input curve */
2752 *out = p->incurve(p, *in, c->iix);
2753 } else if (c->oix >= 0) { /* Output curve */
2754 *out = p->outcurve(p, *in, c->oix);
2755 }
2756 }
2757
2758 /* Function to pass to rspl to set inverse input transfer function, */
2759 /* used for ink limiting calculation. */
2760 static void
icxLuLut_invinput_func(void * cc,double * out,double * in)2761 icxLuLut_invinput_func(
2762 void *cc, /* curvectx structure */
2763 double *out, /* Device output value */
2764 double *in /* Device input value */
2765 ) {
2766 curvectx *c = (curvectx *)cc; /* this */
2767 xfit *p = c->xf;
2768
2769 *out = p->invincurve(p, *in, c->iix);
2770 }
2771
2772
2773 /* Functions to pass to icc settables() to setup icc A2B Lut: */
2774
2775 /* Input table */
set_input(void * cntx,double * out,double * in)2776 static void set_input(void *cntx, double *out, double *in) {
2777 icxLuLut *p = (icxLuLut *)cntx;
2778
2779 if (p->noisluts != 0 && p->noipluts != 0) { /* Input table must be linear */
2780 int i;
2781 for (i = 0; i < p->inputChan; i++)
2782 out[i] = in[i];
2783 } else {
2784 if (p->input(p, out, in) > 1)
2785 error ("%d, %s",p->pp->errc,p->pp->err);
2786 }
2787 }
2788
2789 /* clut */
set_clut(void * cntx,double * out,double * in)2790 static void set_clut(void *cntx, double *out, double *in) {
2791 icxLuLut *p = (icxLuLut *)cntx;
2792 int f;
2793
2794 if (p->clut(p, out, in) > 1)
2795 error ("%d, %s",p->pp->errc,p->pp->err);
2796
2797 /* Convert from efective pcs to natural pcs */
2798 if (p->pcs != p->plu->icp->header->pcs) {
2799 if (p->pcs == icSigLabData)
2800 icmLab2XYZ(&icmD50, out, out);
2801 else
2802 icmXYZ2Lab(&icmD50, out, out);
2803 }
2804 }
2805
2806 /* output */
set_output(void * cntx,double * out,double * in)2807 static void set_output(void *cntx, double *out, double *in) {
2808 icxLuLut *p = (icxLuLut *)cntx;
2809
2810 if (p->nooluts != 0) { /* Output table must be linear */
2811 int i;
2812 for (i = 0; i < p->outputChan; i++)
2813 out[i] = in[i];
2814 } else {
2815 if (p->output(p, out, in) > 1)
2816 error ("%d, %s",p->pp->errc,p->pp->err);
2817 }
2818 }
2819
2820
2821 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2822 /* Routine to figure out a suitable black point for CMYK */
2823
2824 /* Structure to hold optimisation information */
2825 typedef struct {
2826 icxLuLut *p; /* Object being created */
2827 double toAbs[3][3]; /* To abs from aprox relative */
2828 double p1[3]; /* white pivot point in abs Lab */
2829 double p2[3]; /* Point on vector towards black */
2830 double toll; /* Tollerance of black direction */
2831 } bfinds;
2832
2833 /* Optimise device values to minimise L, while remaining */
2834 /* within the ink limit, and staying in line between p1 (white) and p2 (black dir) */
bfindfunc(void * adata,double pv[])2835 static double bfindfunc(void *adata, double pv[]) {
2836 bfinds *b = (bfinds *)adata;
2837 double rv = 0.0;
2838 double tt[3], Lab[3];
2839 co bcc;
2840 double lr, ta, tb, terr; /* L ratio, target a, target b, target error */
2841 double ovr;
2842
2843 //printf("~1 bfindfunc got %f %f %f %f\n", pv[0], pv[1], pv[2], pv[3]);
2844 /* See if over ink limit or outside device range */
2845 ovr = icxLimit(b->p, pv); /* > 0.0 if outside device gamut or ink limit */
2846 if (ovr < 0.0)
2847 ovr = 0.0;
2848 //printf("~1 bfindfunc got ovr %f\n", ovr);
2849
2850 /* Compute the absolute Lab value: */
2851 b->p->input(b->p, bcc.p, pv); /* Through input tables */
2852 b->p->clutTable->interp(b->p->clutTable, &bcc); /* Through clut */
2853 b->p->output(b->p, bcc.v, bcc.v); /* Through the output tables */
2854
2855 if (b->p->pcs != icSigXYZData) /* Convert PCS to XYZ */
2856 icmLab2XYZ(&icmD50, bcc.v, bcc.v);
2857
2858 /* Convert from relative to Absolute colorimetric */
2859 icmMulBy3x3(tt, b->toAbs, bcc.v);
2860 icmXYZ2Lab(&icmD50, Lab, tt); /* Convert to Lab */
2861
2862 #ifdef DEBUG
2863 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]);
2864 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]);
2865 #endif
2866
2867 /* Primary aim is to minimise L value */
2868 rv = Lab[0];
2869
2870 /* See how out of line from p1 to p2 we are */
2871 lr = (Lab[0] - b->p1[0])/(b->p2[0] - b->p1[0]); /* Distance towards p2 from p1 */
2872 ta = lr * (b->p2[1] - b->p1[1]) + b->p1[1]; /* Target a value */
2873 tb = lr * (b->p2[2] - b->p1[2]) + b->p1[2]; /* Target b value */
2874
2875 terr = (ta - Lab[1]) * (ta - Lab[1])
2876 + (tb - Lab[2]) * (tb - Lab[2]);
2877
2878 if (terr < b->toll) /* Tollerance error doesn't count until it's over tollerance */
2879 terr = 0.0;
2880
2881 #ifdef DEBUG
2882 printf("~1 target error %f\n",terr);
2883 #endif
2884 rv += XICC_BLACK_FIND_ABERR_WEIGHT * terr; /* Make ab match more important than min. L */
2885
2886 #ifdef DEBUG
2887 printf("~1 out of range error %f\n",ovr);
2888 #endif
2889 rv += 200 * ovr;
2890
2891 #ifdef DEBUG
2892 printf("~1 black find tc ret %f\n",rv);
2893 #endif
2894 return rv;
2895 }
2896
2897 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2898
2899 /* Create icxLuLut and underlying fwd Lut from scattered data */
2900 /* The scattered data is assumed to map Device -> native PCS */
2901 /* NOTE:- in theory once this icxLuLut is setup, it can be */
2902 /* called to translate color values. In practice I suspect */
2903 /* that the icxLuLut hasn't been setup completely enough to allows this. */
2904 /* Might be easier to close it and re-open it ? */
2905 static icxLuBase *
set_icxLuLut(xicc * xicp,icmLuBase * plu,icmLookupFunc func,icRenderingIntent intent,int flags,int nodp,int nodpbw,cow * ipoints,icxMatrixModel * skm,double dispLuminance,double wpscale,double smooth,double avgdev,double demph,icxViewCond * vc,icxInk * ink,int quality)2906 set_icxLuLut(
2907 xicc *xicp,
2908 icmLuBase *plu, /* Pointer to Lu we are expanding (ours) */
2909 icmLookupFunc func, /* Functionality requested */
2910 icRenderingIntent intent, /* Rendering intent */
2911 int flags, /* white/black point flags etc. */
2912 int nodp, /* Number of points */
2913 int nodpbw, /* Number of points to look for white & black patches in */
2914 cow *ipoints, /* Array of input points (Lab or XYZ normalized to 1.0) */
2915 icxMatrixModel *skm, /* Optional skeleton model (used for input profiles) */
2916 double dispLuminance, /* > 0.0 if display luminance value and is known */
2917 double wpscale, /* > 0.0 if white point is to be scaled */
2918 //double *bpo, /* != NULL for XYZ black point override dev & XYZ */
2919 double smooth, /* RSPL smoothing factor, -ve if raw */
2920 double avgdev, /* reading Average Deviation as a prop. of the input range */
2921 double demph, /* dark emphasis factor for cLUT grid res. */
2922 icxViewCond *vc, /* Viewing Condition (NULL if not using CAM) */
2923 icxInk *ink, /* inking details (NULL for default) */
2924 int quality /* Quality metric, 0..3 */
2925 ) {
2926 icxLuLut *p; /* Object being created */
2927 icc *icco = xicp->pp; /* Underlying icc object */
2928 int luflags = 0; /* icxLuLut alloc clip, merge flags */
2929 int pcsy; /* Effective PCS L or Y chanel index */
2930 double pcsymax; /* Effective PCS L or Y maximum value */
2931 icmHeader *h = icco->header; /* Pointer to icc header */
2932 int maxchan; /* max(inputChan, outputChan) */
2933 int rsplflags = RSPL_NOFLAGS; /* Flags for scattered data rspl */
2934 int e, f, i, j;
2935 double dwhite[MXDI], dblack[MXDI]; /* Device white and black values */
2936 double wp[3]; /* Absolute White point in XYZ */
2937 double bp[3]; /* Absolute Black point in XYZ */
2938 double dgwhite[MXDI]; /* Device space gamut boundary white (ie. RGB 1,1,1) */
2939 double oavgdev[MXDO]; /* Average output value deviation */
2940 int gres[MXDI]; /* RSPL/CLUT resolution */
2941 xfit *xf = NULL; /* Curve fitting class instance */
2942 // co bpop; /* bpo dev + XYZ value */
2943
2944 if (flags & ICX_VERBOSE)
2945 rsplflags |= RSPL_VERBOSE;
2946
2947 // if (flags & ICX_2PASSSMTH)
2948 // rsplflags |= RSPL_2PASSSMTH; /* Smooth data using Gaussian */
2949
2950 // if (flags & ICX_EXTRA_FIT)
2951 // rsplflags |= RSPL_EXTRAFIT2; /* Try extra hard to fit data */
2952
2953 luflags = flags; /* Transfer straight though ? */
2954
2955 /* Do basic creation and initialisation */
2956 if ((p = alloc_icxLuLut(xicp, plu, luflags)) == NULL)
2957 return NULL;
2958
2959 p->func = func;
2960
2961 /* Set LuLut "use" specific creation flags: */
2962 if (flags & ICX_CLIP_NEAREST)
2963 p->nearclip = 1;
2964
2965 /* Set LuLut "create" specific flags: */
2966 if (flags & ICX_NO_IN_SHP_LUTS)
2967 p->noisluts = 1;
2968
2969 if (flags & ICX_NO_IN_POS_LUTS)
2970 p->noipluts = 1;
2971
2972 if (flags & ICX_NO_OUT_LUTS)
2973 p->nooluts = 1;
2974
2975 /* Get the effective spaces of underlying icm, and set icx the same */
2976 plu->spaces(plu, &p->ins, NULL, &p->outs, NULL, NULL, &p->intent, NULL, &p->pcs, NULL);
2977
2978 /* For set_icx the effective pcs has to be the same as the native pcs */
2979
2980 if (p->pcs == icSigXYZData) {
2981 pcsy = 1; /* Y of XYZ */
2982 pcsymax = 1.0;
2983 } else {
2984 pcsy = 0; /* L or Lab */
2985 pcsymax = 100.0;
2986 }
2987
2988 maxchan = p->inputChan > p->outputChan ? p->inputChan : p->outputChan;
2989
2990 /* Translate overall average deviation into output channel deviation */
2991 /* (This is for rspl scattered data fitting smoothness adjustment) */
2992 /* (This could do with more tuning) */
2993
2994 /* XYZ display models are under-smoothed, because the mapping is typically */
2995 /* very "straight", and the lack of tension reduces any noise reduction effect. */
2996 /* !!! This probably means that we should switch to 3rd order smoothness criteria !! */
2997 /* We apply an arbitrary correction here */
2998 /* !!!! There is also a bug in the rspl code, where smoothness is */
2999 /* scaled by data range. This is making Lab smoothing ~100 times */
3000 /* more than XYZ smoothing. Fix this with SMOOTH2 changes ?? */
3001 if (p->pcs == icSigXYZData) {
3002 oavgdev[0] = XYZ_EXTRA_SMOOTH * 0.70 * avgdev;
3003 oavgdev[1] = XYZ_EXTRA_SMOOTH * 1.00 * avgdev;
3004 oavgdev[2] = XYZ_EXTRA_SMOOTH * 0.70 * avgdev;
3005 } else if (p->pcs == icSigLabData) {
3006 oavgdev[0] = 1.00 * avgdev;
3007 oavgdev[1] = 0.70 * avgdev;
3008 oavgdev[2] = 0.70 * avgdev;
3009 } else { /* Hmmm */
3010 for (f = 0; f < p->outputChan; f++)
3011 oavgdev[f] = avgdev;
3012 }
3013
3014 /* In general the native and effective ranges of the icx will be the same as the */
3015 /* underlying icm lookup object. */
3016 p->plu->get_lutranges(p->plu, p->ninmin, p->ninmax, p->noutmin, p->noutmax);
3017 p->plu->get_ranges(p->plu, p->inmin, p->inmax, p->outmin, p->outmax);
3018
3019 /* ??? Does this ever happen with set_icxLuLut() ??? */
3020 /* If we have a Jab PCS override, reflect this in the effective icx range. */
3021 /* Note that the ab ranges are nominal. They will exceed this range */
3022 /* for colors representable in L*a*b* PCS */
3023 if (p->ins == icxSigJabData) {
3024 p->inmin[0] = 0.0; p->inmax[0] = 100.0;
3025 p->inmin[1] = -128.0; p->inmax[1] = 128.0;
3026 p->inmin[2] = -128.0; p->inmax[2] = 128.0;
3027 } else if (p->outs == icxSigJabData) {
3028 p->outmin[0] = 0.0; p->outmax[0] = 100.0;
3029 p->outmin[1] = -128.0; p->outmax[1] = 128.0;
3030 p->outmin[2] = -128.0; p->outmax[2] = 128.0;
3031 }
3032
3033 /* ------------------------------- */
3034
3035 if (flags & ICX_VERBOSE)
3036 printf("Estimating white point\n");
3037
3038 icmXYZ2Ary(wp, icmD50); /* Set a default value - D50 */
3039 icmXYZ2Ary(bp, icmBlack); /* Set a default value - absolute black */
3040
3041 {
3042 /* Figure out as best we can the device white and black points */
3043
3044 /* Compute device white and black points as if */
3045 /* we are doing an Output or Display device */
3046 {
3047 switch (h->colorSpace) {
3048
3049 case icSigCmykData:
3050 for (e = 0; e < p->inputChan; e++) {
3051 dwhite[e] = 0.0;
3052 dblack[e] = 1.0;
3053 }
3054 break;
3055 case icSigCmyData:
3056 for (e = 0; e < p->inputChan; e++) {
3057 dwhite[e] = 0.0;
3058 dblack[e] = 1.0;
3059 }
3060 break;
3061 case icSigRgbData:
3062 for (e = 0; e < p->inputChan; e++) {
3063 dwhite[e] = 1.0;
3064 dblack[e] = 0.0;
3065 }
3066 break;
3067
3068 case icSigGrayData: { /* Could be additive or subtractive */
3069 double maxY, minY; /* Y min and max */
3070 double maxd, mind; /* Corresponding device values */
3071
3072 maxY = -1e8, minY = 1e8;
3073
3074 /* Figure out if it's additive or subtractive */
3075 for (i = 0; i < nodpbw; i++) {
3076 double Y;
3077 if (p->pcs != icSigXYZData) { /* Convert white point to XYZ */
3078 double xyz[3];
3079 icmLab2XYZ(&icmD50, xyz, ipoints[i].v);
3080 Y = xyz[1];
3081 } else {
3082 Y = ipoints[i].v[1];
3083 }
3084
3085 if (Y > maxY) {
3086 maxY = Y;
3087 maxd = ipoints[i].p[0];
3088 }
3089 if (Y < minY) {
3090 minY = Y;
3091 mind = ipoints[i].p[0];
3092 }
3093 }
3094 if (maxd < mind) { /* Subtractive */
3095 for (e = 0; e < p->inputChan; e++) {
3096 dwhite[e] = 0.0;
3097 dblack[e] = 1.0;
3098 }
3099 } else { /* Additive */
3100 for (e = 0; e < p->inputChan; e++) {
3101 dwhite[e] = 1.0;
3102 dblack[e] = 0.0;
3103 }
3104 }
3105 break;
3106 }
3107 default:
3108 xicp->errc = 1;
3109 sprintf(xicp->err,"set_icxLuLut: can't handle color space %s",
3110 icm2str(icmColorSpaceSignature, h->colorSpace));
3111 p->del((icxLuBase *)p);
3112 return NULL;
3113 break;
3114 }
3115 }
3116 /* dwhite is what we want for dgwhite[], used for XFIT_OUT_WP_REL_US */
3117 for (e = 0; e < p->inputChan; e++)
3118 dgwhite[e] = dwhite[e];
3119
3120 /* If this is actuall an input device, lookup wp & bp */
3121 /* and override dwhite & dblack */
3122 if (h->deviceClass == icSigInputClass) {
3123 double wpy = -1e60, bpy = 1e60;
3124 int wix = -1, bix = -1;
3125 /* We assume that the input target is well behaved, */
3126 /* and that it includes a white and black point patch, */
3127 /* and that they have the extreme L/Y values */
3128
3129 /*
3130 NOTE that this may not be the best approach !
3131 It may be better to average the chromaticity
3132 of all the neutral seeming patches, since
3133 the whitest patch may have (for instance)
3134 a blue tint.
3135 */
3136
3137 /* Discover the white and black patches */
3138 for (i = 0; i < nodpbw; i++) {
3139 double labv[3], yv;
3140
3141 /* Create D50 Lab to allow some chromatic sensitivity */
3142 /* in picking the white point */
3143 if (p->pcs == icSigXYZData)
3144 icmXYZ2Lab(&icmD50, labv, ipoints[i].v);
3145 else
3146 icmCpy3(labv,ipoints[i].v);
3147
3148 #ifdef NEVER
3149 /* Choose largest Y or L* */
3150 if (ipoints[i].v[pcsy] > wpy) {
3151 wp[0] = ipoints[i].v[0];
3152 wp[1] = ipoints[i].v[1];
3153 wp[2] = ipoints[i].v[2];
3154 for (e = 0; e < p->inputChan; e++)
3155 dwhite[e] = ipoints[i].p[e];
3156 wpy = ipoints[i].v[pcsy];
3157 wix = i;
3158 }
3159 #else
3160
3161 /* Tilt things towards D50 neutral white patches */
3162 yv = labv[0] - 0.3 * sqrt(labv[1] * labv[1] + labv[2] * labv[2]);
3163 //printf("~1 patch %d Lab = %s, yv = %f\n",i+1,icmPdv(3,labv),yv);
3164 if (yv > wpy) {
3165 //printf("~1 best so far\n");
3166 wp[0] = ipoints[i].v[0];
3167 wp[1] = ipoints[i].v[1];
3168 wp[2] = ipoints[i].v[2];
3169 for (e = 0; e < p->inputChan; e++)
3170 dwhite[e] = ipoints[i].p[e];
3171 wpy = yv;
3172 wix = i;
3173 }
3174 #endif
3175 if (ipoints[i].v[pcsy] < bpy) {
3176 bp[0] = ipoints[i].v[0];
3177 bp[1] = ipoints[i].v[1];
3178 bp[2] = ipoints[i].v[2];
3179 for (e = 0; e < p->inputChan; e++)
3180 dblack[e] = ipoints[i].p[e];
3181 bpy = ipoints[i].v[pcsy];
3182 bix = i;
3183 }
3184 }
3185 /* Make sure values are XYZ */
3186 if (p->pcs != icSigXYZData) {
3187 icmLab2XYZ(&icmD50, wp, wp);
3188 icmLab2XYZ(&icmD50, bp, bp);
3189 }
3190
3191 if (flags & ICX_VERBOSE) {
3192 printf("Picked white patch %d with dev = %s\n XYZ = %s, Lab = %s\n",
3193 wix+1, icmPdv(p->inputChan, dwhite), icmPdv(3, wp), icmPLab(wp));
3194 printf("Picked black patch %d with dev = %s\n XYZ = %s, Lab = %s\n",
3195 bix+1, icmPdv(p->inputChan, dblack), icmPdv(3, bp), icmPLab(bp));
3196 }
3197
3198 /* else Output or Display device */
3199 } else {
3200 /* We assume that the output target is well behaved, */
3201 /* and that it includes a white point patch. */
3202 int nw = 0;
3203
3204 wp[0] = wp[1] = wp[2] = 0.0;
3205
3206 switch (h->colorSpace) {
3207
3208 case icSigCmykData:
3209 for (i = 0; i < nodpbw; i++) {
3210 if (ipoints[i].p[0] < 0.001
3211 && ipoints[i].p[1] < 0.001
3212 && ipoints[i].p[2] < 0.001
3213 && ipoints[i].p[3] < 0.001) {
3214 wp[0] += ipoints[i].v[0];
3215 wp[1] += ipoints[i].v[1];
3216 wp[2] += ipoints[i].v[2];
3217 nw++;
3218 }
3219 }
3220 break;
3221 case icSigCmyData:
3222 for (i = 0; i < nodpbw; i++) {
3223 if (ipoints[i].p[0] < 0.001
3224 && ipoints[i].p[1] < 0.001
3225 && ipoints[i].p[2] < 0.001) {
3226 wp[0] += ipoints[i].v[0];
3227 wp[1] += ipoints[i].v[1];
3228 wp[2] += ipoints[i].v[2];
3229 nw++;
3230 }
3231 }
3232 break;
3233 case icSigRgbData:
3234 for (i = 0; i < nodpbw; i++) {
3235 if (ipoints[i].p[0] > 0.999
3236 && ipoints[i].p[1] > 0.999
3237 && ipoints[i].p[2] > 0.999) {
3238 wp[0] += ipoints[i].v[0];
3239 wp[1] += ipoints[i].v[1];
3240 wp[2] += ipoints[i].v[2];
3241 nw++;
3242 }
3243 }
3244 /* Setup bpo device value in case we need it */
3245 // bpop.p[0] = bpop.p[1] = bpop.p[2] = 0.0;
3246 break;
3247
3248 case icSigGrayData: { /* Could be additive or subtractive */
3249 double minwp[3], maxwp[3];
3250 int nminwp = 0, nmaxwp = 0;
3251
3252 minwp[0] = minwp[1] = minwp[2] = 0.0;
3253 maxwp[0] = maxwp[1] = maxwp[2] = 0.0;
3254
3255 /* Look for both */
3256 for (i = 0; i < nodpbw; i++) {
3257 if (ipoints[i].p[0] < 0.001)
3258 minwp[0] += ipoints[i].v[0];
3259 minwp[1] += ipoints[i].v[1];
3260 minwp[2] += ipoints[i].v[2]; {
3261 nminwp++;
3262 }
3263 if (ipoints[i].p[0] > 0.999)
3264 maxwp[0] += ipoints[i].v[0];
3265 maxwp[1] += ipoints[i].v[1];
3266 maxwp[2] += ipoints[i].v[2]; {
3267 nmaxwp++;
3268 }
3269 }
3270 if (nminwp > 0) { /* Subtractive */
3271 wp[0] = minwp[0];
3272 wp[1] = minwp[1];
3273 wp[2] = minwp[2];
3274 nw = nminwp;
3275 if (minwp[pcsy]/nminwp < (0.5 * pcsymax))
3276 nw = 0; /* Looks like a mistake */
3277 // bpop.p[0] = 1.0;
3278 }
3279 if (nmaxwp > 0 /* Additive */
3280 && (nminwp == 0 || maxwp[pcsy]/nmaxwp > minwp[pcsy]/nminwp)) {
3281 wp[0] = maxwp[0];
3282 wp[1] = maxwp[1];
3283 wp[2] = maxwp[2];
3284 nw = nmaxwp;
3285 if (maxwp[pcsy]/nmaxwp < (0.5 * pcsymax))
3286 nw = 0; /* Looks like a mistake */
3287 // bpop.p[0] = 0.0;
3288 }
3289 break;
3290 }
3291
3292 default:
3293 xicp->errc = 1;
3294 sprintf(xicp->err,"set_icxLuLut: can't handle color space %s",
3295 icm2str(icmColorSpaceSignature, h->colorSpace));
3296 p->del((icxLuBase *)p);
3297 return NULL;
3298 break;
3299 }
3300
3301 if (nw == 0) {
3302 xicp->errc = 1;
3303 sprintf(xicp->err,"set_icxLuLut: can't handle test points without a white patch");
3304 p->del((icxLuBase *)p);
3305 return NULL;
3306 }
3307 wp[0] /= (double)nw;
3308 wp[1] /= (double)nw;
3309 wp[2] /= (double)nw;
3310 if (p->pcs != icSigXYZData) /* Convert white point to XYZ */
3311 icmLab2XYZ(&icmD50, wp, wp);
3312
3313 // if (bpo != NULL) { /* Copy black override XYZ value */
3314 // bpop.v[0] = bpo[0];
3315 // bpop.v[1] = bpo[1];
3316 // bpop.v[2] = bpo[2];
3317 // }
3318 }
3319
3320 if (flags & ICX_VERBOSE) {
3321 printf("Approximate White point XYZ = %s, Lab = %s\n", icmPdv(3,wp),icmPLab(wp));
3322 }
3323 }
3324
3325 if (h->colorSpace == icSigGrayData) { /* Don't use device or PCS curves for monochrome */
3326 p->noisluts = p->noipluts = p->nooluts = 1;
3327 }
3328
3329 if ((flags & ICX_VERBOSE) && (p->noisluts == 0 || p->noipluts == 0 || p->nooluts == 0))
3330 printf("Creating optimised per channel curves\n");
3331
3332 /* Set the target CLUT grid resolution so in/out curves can be optimised for it */
3333 for (e = 0; e < p->inputChan; e++)
3334 gres[e] = p->lut->clutPoints;
3335
3336 /* Setup and then create xfit object that does most of the work */
3337 {
3338 int xfflags = 0; /* xfit flags */
3339 double in_min[MXDI]; /* Input value scaling minimum */
3340 double in_max[MXDI]; /* Input value scaling maximum */
3341 double out_min[MXDO]; /* Output value scaling minimum */
3342 double out_max[MXDO]; /* Output value scaling maximum */
3343 int iluord, sluord, oluord;
3344 int iord[MXDI]; /* Input curve orders */
3345 int sord[MXDI]; /* Input sub-grid curve orders */
3346 int oord[MXDO]; /* Output curve orders */
3347 double shp_smooth[MXDI];/* Smoothing factors for each curve, nom = 1.0 */
3348 double out_smooth[MXDO];
3349
3350 optcomb tcomb = oc_ipo; /* Create all by default */
3351
3352 if ((xf = new_xfit(icco)) == NULL) {
3353 p->pp->errc = 2;
3354 sprintf(p->pp->err,"Creation of xfit object failed");
3355 p->del((icxLuBase *)p);
3356 return NULL;
3357 }
3358
3359 /* Setup for optimising run */
3360 if (p->noisluts)
3361 tcomb &= ~oc_i;
3362
3363 if (p->noipluts)
3364 tcomb &= ~oc_p;
3365
3366 if (p->nooluts)
3367 tcomb &= ~oc_o;
3368
3369 if (flags & ICX_VERBOSE)
3370 xfflags |= XFIT_VERB;
3371
3372 if (flags & ICX_SET_WHITE) {
3373
3374 xfflags |= XFIT_OUT_WP_REL;
3375 if ((flags & ICX_SET_WHITE_C) == ICX_SET_WHITE_C) {
3376 xfflags |= XFIT_OUT_WP_REL_C;
3377 }
3378 else if ((flags & ICX_SET_WHITE_US) == ICX_SET_WHITE_US)
3379 xfflags |= XFIT_OUT_WP_REL_US;
3380
3381 if (p->pcs != icSigXYZData)
3382 xfflags |= XFIT_OUT_LAB;
3383 }
3384
3385 /* If asked, clip the absolute white point to have Y <= 1.0 ? */
3386 if (flags & ICX_CLIP_WB)
3387 xfflags |= XFIT_CLIP_WP;
3388
3389 /* With current B2A code, make sure a & b curves */
3390 /* pass through zero. */
3391 if (p->pcs == icSigLabData) {
3392 xfflags |= XFIT_OUT_ZERO;
3393 }
3394
3395 /* Let xfit create the clut */
3396 xfflags |= XFIT_MAKE_CLUT;
3397
3398 /* Set the curve orders for input (device) and output (PCS) */
3399 if (quality >= 3) { /* Ultra high */
3400 iluord = 25;
3401 sluord = 4;
3402 oluord = 25;
3403 } else if (quality == 2) { /* High */
3404 iluord = 20;
3405 sluord = 3;
3406 oluord = 20;
3407 } else if (quality == 1) { /* Medium */
3408 iluord = 17;
3409 sluord = 2;
3410 oluord = 17;
3411 } else { /* Low */
3412 iluord = 10;
3413 sluord = 1;
3414 oluord = 10;
3415 }
3416 for (e = 0; e < p->inputChan; e++) {
3417 iord[e] = iluord;
3418 sord[e] = sluord;
3419 in_min[e] = p->inmin[e];
3420 in_max[e] = p->inmax[e];
3421 shp_smooth[e] = SHP_SMOOTH;
3422 }
3423
3424 for (f = 0; f < p->outputChan; f++) {
3425 oord[f] = oluord;
3426 out_min[f] = p->outmin[f];
3427 out_max[f] = p->outmax[f];
3428
3429 /* Hack to prevent a convex L curve pushing */
3430 /* the clut L values above the maximum value */
3431 /* that can be represented, causing clipping. */
3432 /* Do this by making sure that the L curve pivots */
3433 /* through 100.0 to 100.0 */
3434 if (f == 0 && p->pcs == icSigLabData) {
3435 if (out_min[f] < 0.0001 && out_max[f] > 100.0) {
3436 out_max[f] = 100.0;
3437 }
3438 }
3439 out_smooth[f] = OUT_SMOOTH1;
3440 if (f != 0 && p->pcs == icSigLabData) /* a* & b* */
3441 out_smooth[f] = OUT_SMOOTH2;
3442 }
3443
3444 //printf("~1 wp before xfit %f %f %f\n", wp[0], wp[1], wp[2]);
3445 /* Create input, sub and output per channel curves (if configured), */
3446 /* adjust for white point to make relative (if configured), */
3447 /* and create clut rspl using xfit class. */
3448 /* The true white point for the returned curves and rspl is returned. */
3449 if (xf->fit(xf, xfflags, p->inputChan, p->outputChan,
3450 rsplflags, wp, dwhite, wpscale, dgwhite,
3451 ipoints, nodp, skm, in_min, in_max, gres, out_min, out_max,
3452 // bpo != NULL ? &bpop : NULL,
3453 smooth, oavgdev, demph, iord, sord, oord, shp_smooth, out_smooth, tcomb,
3454 (void *)p, xfit_to_de2, xfit_to_dde2) != 0) {
3455 p->pp->errc = 2;
3456 sprintf(p->pp->err,"xfit fitting failed");
3457 xf->del(xf);
3458 p->del((icxLuBase *)p);
3459 return NULL;
3460 }
3461 //printf("~1 wp after xfit %f %f %f\n", wp[0], wp[1], wp[2]);
3462
3463 /* - - - - - - - - - - - - - - - */
3464 /* Set the xicc input curve rspl */
3465 for (e = 0; e < p->inputChan; e++) {
3466 curvectx cx;
3467
3468 cx.xf = xf;
3469 cx.oix = -1;
3470 cx.iix = e;
3471
3472 if ((p->inputTable[e] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) {
3473 p->pp->errc = 2;
3474 sprintf(p->pp->err,"Creation of input table rspl failed");
3475 xf->del(xf);
3476 p->del((icxLuBase *)p);
3477 return NULL;
3478 }
3479
3480 p->inputTable[e]->set_rspl(p->inputTable[e], RSPL_NOFLAGS,
3481 (void *)&cx, set_linfunc,
3482 &p->ninmin[e], &p->ninmax[e],
3483 (int *)&p->lut->inputEnt,
3484 &p->ninmin[e], &p->ninmax[e]);
3485 }
3486
3487 /* - - - - - - - - - - - - - - - */
3488 /* Set the xicc output curve rspl */
3489 for (f = 0; f < p->outputChan; f++) {
3490 curvectx cx;
3491
3492 cx.xf = xf;
3493 cx.iix = -1;
3494 cx.oix = f;
3495
3496 if ((p->outputTable[f] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) {
3497 p->pp->errc = 2;
3498 sprintf(p->pp->err,"Creation of output table rspl failed");
3499 xf->del(xf);
3500 p->del((icxLuBase *)p);
3501 return NULL;
3502 }
3503
3504 p->outputTable[f]->set_rspl(p->outputTable[f], RSPL_NOFLAGS,
3505 (void *)&cx, set_linfunc,
3506 &p->noutmin[f], &p->noutmax[f],
3507 (int *)&p->lut->outputEnt,
3508 &p->noutmin[f], &p->noutmax[f]);
3509
3510 }
3511 }
3512
3513 if (flags & ICX_VERBOSE)
3514 printf("Creating fast inverse input lookups\n");
3515
3516 /* Create rspl based reverse input lookups used in ink limit function. */
3517 for (e = 0; e < p->inputChan; e++) {
3518 int res = 256;
3519 curvectx cx;
3520
3521 cx.xf = xf;
3522 cx.oix = -1;
3523 cx.iix = e;
3524
3525 if ((p->revinputTable[e] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) {
3526 p->pp->errc = 2;
3527 sprintf(p->pp->err,"Creation of reverse input table rspl failed");
3528 xf->del(xf);
3529 p->del((icxLuBase *)p);
3530 return NULL;
3531 }
3532 p->revinputTable[e]->set_rspl(p->revinputTable[e], RSPL_NOFLAGS,
3533 (void *)&cx, icxLuLut_invinput_func,
3534 &p->ninmin[e], &p->ninmax[e], &res, &p->ninmin[e], &p->ninmax[e]);
3535 }
3536
3537
3538 /* ------------------------------- */
3539 /* Set clut lookup table from xfit */
3540 p->clutTable = xf->clut;
3541 xf->clut = NULL;
3542
3543 /* Setup all the clipping, ink limiting and auxiliary stuff, */
3544 /* in case a reverse call is used. Need to avoid relying on inking */
3545 /* stuff that makes use of the white/black points, since they haven't */
3546 /* been set up properly yet. */
3547 if (setup_ink_icxLuLut(p, ink, 0) != 0) {
3548 xf->del(xf);
3549 p->del((icxLuBase *)p);
3550 return NULL;
3551 }
3552
3553 /* Deal with finalizing white/black points */
3554 {
3555
3556 if ((flags & ICX_SET_WHITE) && (flags & ICX_VERBOSE)) {
3557 printf("White point XYZ = %s, Lab = %s\n", icmPdv(3,wp),icmPLab(wp));
3558 }
3559
3560 /* Lookup the black point */
3561 { /* Black Point Tag: */
3562 co bcc;
3563
3564 if (flags & ICX_VERBOSE)
3565 printf("Find black point\n");
3566
3567 /* !!! Hmm. For CMY and RGB we are simply using the device */
3568 /* combination values as the black point. In reality we might */
3569 /* want to have the option of using a neutral black point, */
3570 /* just like CMYK ?? */
3571
3572 /* For CMYK devices, we choose a black that has minumum L within */
3573 /* the ink limits, and if XICC_NEUTRAL_CMYK_BLACK it will in the direction */
3574 /* that has the same chromaticity as the white point, else choose the same */
3575 /* Lab vector direction as K, with the minimum L value. */
3576 /* (Note this is duplicated in xicc.c icxLu_comp_bk_point() !!!) */
3577 if (h->deviceClass != icSigInputClass
3578 && h->colorSpace == icSigCmykData) {
3579 bfinds bfs; /* Callback context */
3580 double sr[MXDO]; /* search radius */
3581 double tt[MXDO]; /* Temporary */
3582 double rs0[MXDO], rs1[MXDO]; /* Random start candidates */
3583 int trial;
3584 double brv;
3585 int kch = 3;
3586
3587 /* Setup callback function context */
3588 bfs.p = p;
3589 bfs.toll = XICC_BLACK_POINT_TOLL;
3590
3591 /* !!! we should use an accessor funcion of xfit !!! */
3592 for (i = 0; i < 3; i++) {
3593 for (j = 0; j < 3; j++) {
3594 bfs.toAbs[i][j] = xf->toAbs[i][j];
3595 }
3596 }
3597
3598 /* Lookup abs Lab value of white point */
3599 icmXYZ2Lab(&icmD50, bfs.p1, wp);
3600
3601 #ifdef XICC_NEUTRAL_CMYK_BLACK
3602 icmScale3(tt, wp, 0.02); /* Scale white XYZ towards 0,0,0 */
3603 icmXYZ2Lab(&icmD50, bfs.p2, tt); /* Convert black direction to Lab */
3604 if (flags & ICX_VERBOSE)
3605 printf("Neutral black direction (Lab) = %f %f %f\n",bfs.p2[0], bfs.p2[1], bfs.p2[2]);
3606
3607 #else /* K direction black */
3608 /* Now figure abs Lab value of K only, as the direction */
3609 /* to use for the rich black. */
3610 for (e = 0; e < p->inputChan; e++)
3611 bcc.p[e] = 0.0;
3612 if (p->ink.klimit <= 0.0)
3613 bcc.p[kch] = 1.0;
3614 else
3615 bcc.p[kch] = p->ink.klimit; /* K value */
3616
3617 p->input(p, bcc.p, bcc.p); /* Through input tables */
3618 p->clutTable->interp(p->clutTable, &bcc); /* Through clut */
3619 p->output(p, bcc.v, bcc.v); /* Through the output tables */
3620
3621 if (p->pcs != icSigXYZData) /* Convert PCS to XYZ */
3622 icmLab2XYZ(&icmD50, bcc.v, bcc.v);
3623
3624 /* Convert from relative to Absolute colorimetric */
3625 icmMulBy3x3(tt, xf->toAbs, bcc.v);
3626 icmXYZ2Lab(&icmD50, bfs.p2, tt); /* Convert K only black point to Lab */
3627
3628 if (flags & ICX_VERBOSE)
3629 printf("K only black direction (Lab) = %f %f %f\n",bfs.p2[0], bfs.p2[1], bfs.p2[2]);
3630 #endif
3631 /* Set the random start 0 location as 000K */
3632 /* and the random start 1 location as CMY0 */
3633 {
3634 double tv;
3635
3636 for (e = 0; e < p->inputChan; e++)
3637 rs0[e] = 0.0;
3638 if (p->ink.klimit <= 0.0)
3639 rs0[kch] = 1.0;
3640 else
3641 rs0[kch] = p->ink.klimit; /* K value */
3642
3643 if (p->ink.tlimit <= 0.0)
3644 tv = 1.0;
3645 else
3646 tv = p->ink.tlimit/(p->inputChan - 1.0);
3647 for (e = 0; e < p->inputChan; e++)
3648 rs1[e] = tv;
3649 rs1[kch] = 0.0; /* K value */
3650 }
3651
3652 /* Start with the K only as the current best value */
3653 for (e = 0; e < p->inputChan; e++)
3654 bcc.p[e] = rs0[e];
3655 brv = bfindfunc((void *)&bfs, bcc.p);
3656 //printf("~1 initial brv for K only = %f\n",brv);
3657
3658 /* Find the device black point using optimization */
3659 /* Do several trials to avoid local minima. */
3660 rand32(0x12345678); /* Make trial values deterministic */
3661 for (trial = 0; trial < 200; trial++) {
3662 double rv; /* Temporary */
3663
3664 /* Start first trial at 000K */
3665 if (trial == 0) {
3666 for (j = 0; j < p->inputChan; j++) {
3667 tt[j] = rs0[j];
3668 sr[j] = 0.1;
3669 }
3670
3671 } else {
3672 /* Base is random between 000K and CMY0: */
3673 if (trial < 100) {
3674 rv = d_rand(0.0, 1.0);
3675 for (j = 0; j < p->inputChan; j++) {
3676 tt[j] = rv * rs0[j] + (1.0 - rv) * rs1[j];
3677 sr[j] = 0.1;
3678 }
3679 /* Base on current best */
3680 } else {
3681 for (j = 0; j < p->inputChan; j++) {
3682 tt[j] = bcc.p[j];
3683 sr[j] = 0.1;
3684 }
3685 }
3686
3687 /* Then add random start offset */
3688 for (rv = 0.0, j = 0; j < p->inputChan; j++) {
3689 tt[j] += d_rand(-0.5, 0.5);
3690 if (tt[j] < 0.0)
3691 tt[j] = 0.0;
3692 else if (tt[j] > 1.0)
3693 tt[j] = 1.0;
3694 }
3695 }
3696
3697 /* Clip to be within ink limit */
3698 icxDoLimit(p, tt, tt);
3699
3700 if (powell(&rv, p->inputChan, tt, sr, 0.000001, 1000,
3701 bfindfunc, (void *)&bfs, NULL, NULL) == 0) {
3702 //printf("~1 trial %d, rv %f bp %f %f %f %f\n",trial,rv,tt[0],tt[1],tt[2],tt[3]);
3703 if (rv < brv) {
3704 //printf("~1 new best\n");
3705 brv = rv;
3706 for (j = 0; j < p->inputChan; j++)
3707 bcc.p[j] = tt[j];
3708 }
3709 }
3710 }
3711 if (brv > 1000.0)
3712 error ("set_icxLuLut: Black point powell failed");
3713
3714 /* Make sure resulting device values are strictly in range */
3715 for (j = 0; j < p->inputChan; j++) {
3716 if (bcc.p[j] < 0.0)
3717 bcc.p[j] = 0.0;
3718 else if (bcc.p[j] > 1.0)
3719 bcc.p[j] = 1.0;
3720 }
3721 /* Now have device black in bcc.p[] */
3722 //printf("~1 Black point is CMYK %f %f %f %f\n", bcc.p[0], bcc.p[1], bcc.p[2], bcc.p[3]);
3723
3724 /* Else not a CMYK output device, */
3725 /* use the previously determined device black value. */
3726 } else {
3727 for (e = 0; e < p->inputChan; e++)
3728 bcc.p[e] = dblack[e];
3729 }
3730
3731 /* Lookup the PCS for the device black: */
3732 p->input(p, bcc.p, bcc.p); /* Through input tables */
3733 p->clutTable->interp(p->clutTable, &bcc); /* Through clut */
3734 p->output(p, bcc.v, bcc.v); /* Through the output tables */
3735
3736 if (p->pcs != icSigXYZData) /* Convert PCS to XYZ */
3737 icmLab2XYZ(&icmD50, bcc.v, bcc.v);
3738
3739 /* Convert from relative to Absolute colorimetric */
3740 icmMulBy3x3(bp, xf->toAbs, bcc.v);
3741
3742 /* Got XYZ black point in bp[] */
3743 if (flags & ICX_VERBOSE) {
3744 printf("Black point XYZ = %s, Lab = %s\n", icmPdv(3,bp),icmPLab(bp));
3745 }
3746
3747 /* Some ICC sw gets upset if the bp is at all -ve. */
3748 /* Clip it if this is the case */
3749 /* (Or we could get xfit to rescale the rspl instead ??) */
3750 if ((flags & ICX_CLIP_WB)
3751 && (bp[0] < 0.0 || bp[1] < 0.0 || bp[2] < 0.0)) {
3752 if (bp[0] < 0.0)
3753 bp[0] = 0.0;
3754 if (bp[1] < 0.0)
3755 bp[1] = 0.0;
3756 if (bp[2] < 0.0)
3757 bp[2] = 0.0;
3758 if (flags & ICX_VERBOSE) {
3759 printf("Black point clipped to XYZ = %s, Lab = %s\n",icmPdv(3,bp),icmPLab(bp));
3760 }
3761 }
3762 }
3763
3764 /* If this is a display, adjust the white point to be */
3765 /* exactly Y = 1.0, and compensate the dispLuminance */
3766 /* and black point accordingly. The Lut is already set to */
3767 /* assume device white maps to perfect PCS white. */
3768 if (h->deviceClass == icSigDisplayClass) {
3769 double scale = 1.0/wp[1];
3770 int i;
3771
3772 dispLuminance /= scale;
3773
3774 for (i = 0; i < 3; i++) {
3775 wp[i] *= scale;
3776 bp[i] *= scale;
3777 }
3778
3779 // if (bpo != NULL) {
3780 // bp[0] = bpo[0];
3781 // bp[1] = bpo[1];
3782 // bp[2] = bpo[2];
3783 // printf("Overide Black point XYZ = %s, Lab = %s\n", icmPdv(3,bp),icmPLab(bp));
3784 // }
3785 }
3786
3787 if (h->deviceClass == icSigDisplayClass
3788 && dispLuminance > 0.0) {
3789 icmXYZArray *wo;
3790 if ((wo = (icmXYZArray *)icco->read_tag(
3791 icco, icSigLuminanceTag)) == NULL) {
3792 xicp->errc = 1;
3793 sprintf(xicp->err,"icx_set_luminance: couldn't find luminance tag");
3794 p->del((icxLuBase *)p);
3795 return NULL;
3796 }
3797 if (wo->ttype != icSigXYZArrayType) {
3798 xicp->errc = 1;
3799 sprintf(xicp->err,"luminance: tag has wrong type");
3800 p->del((icxLuBase *)p);
3801 return NULL;
3802 }
3803
3804 wo->size = 1;
3805 wo->allocate((icmBase *)wo); /* Allocate space */
3806 wo->data[0].X = 0.0;
3807 wo->data[0].Y = dispLuminance * wp[1];
3808 wo->data[0].Z = 0.0;
3809
3810 if (flags & ICX_VERBOSE)
3811 printf("Display Luminance = %f\n", wo->data[0].Y);
3812 }
3813
3814 /* Write white and black points */
3815 if (flags & ICX_SET_WHITE) { /* White Point Tag: */
3816 icmXYZArray *wo;
3817 if ((wo = (icmXYZArray *)icco->read_tag(
3818 icco, icSigMediaWhitePointTag)) == NULL) {
3819 xicp->errc = 1;
3820 sprintf(xicp->err,"icx_set_white_black: couldn't find white tag");
3821 xf->del(xf);
3822 p->del((icxLuBase *)p);
3823 return NULL;
3824 }
3825 if (wo->ttype != icSigXYZArrayType) {
3826 xicp->errc = 1;
3827 sprintf(xicp->err,"icx_set_white_black: white tag has wrong type");
3828 xf->del(xf);
3829 p->del((icxLuBase *)p);
3830 return NULL;
3831 }
3832
3833 wo->size = 1;
3834 wo->allocate((icmBase *)wo); /* Allocate space */
3835 wo->data[0].X = wp[0];
3836 wo->data[0].Y = wp[1];
3837 wo->data[0].Z = wp[2];
3838 }
3839 if (flags & ICX_SET_BLACK) { /* Black Point Tag: */
3840 icmXYZArray *wo;
3841 if ((wo = (icmXYZArray *)icco->read_tag(
3842 icco, icSigMediaBlackPointTag)) == NULL) {
3843 xicp->errc = 1;
3844 sprintf(xicp->err,"icx_set_white_black: couldn't find black tag");
3845 xf->del(xf);
3846 p->del((icxLuBase *)p);
3847 return NULL;
3848 }
3849 if (wo->ttype != icSigXYZArrayType) {
3850 xicp->errc = 1;
3851 sprintf(xicp->err,"icx_set_white_black: black tag has wrong type");
3852 xf->del(xf);
3853 p->del((icxLuBase *)p);
3854 return NULL;
3855 }
3856
3857 wo->size = 1;
3858 wo->allocate((icmBase *)wo); /* Allocate space */
3859 wo->data[0].X = bp[0];
3860 wo->data[0].Y = bp[1];
3861 wo->data[0].Z = bp[2];
3862 }
3863 if ((flags & ICX_SET_WHITE) || (flags & ICX_SET_BLACK)) {
3864 /* Make sure ICC white/black point lookup notices the new white and black points */
3865 p->plu->init_wh_bk(p->plu);
3866 }
3867
3868 /* Setup the clut clipping, ink limiting and auxiliary stuff again */
3869 /* since re_set_rspl will have invalidated */
3870 if (setup_ink_icxLuLut(p, ink, 1) != 0) {
3871 xf->del(xf);
3872 p->del((icxLuBase *)p);
3873 return NULL;
3874 }
3875 }
3876
3877 /* Done with xfit now */
3878 xf->del(xf);
3879 xf = NULL;
3880
3881 if (setup_clip_icxLuLut(p) != 0) {
3882 p->del((icxLuBase *)p);
3883 return NULL;
3884 }
3885
3886 /* ------------------------------- */
3887
3888 //Debugging clipping of clut
3889 //printf("~1 xlut.c: noutmin = %f %f %f\n", p->noutmin[0], p->noutmin[1], p->noutmin[2]);
3890 //printf("~1 xlut.c: noutmax = %f %f %f\n", p->noutmax[0], p->noutmax[1], p->noutmax[2]);
3891 //printf("~1 xlut.c: outmin = %f %f %f\n", p->outmin[0], p->outmin[1], p->outmin[2]);
3892 //printf("~1 xlut.c: outmax = %f %f %f\n", p->outmax[0], p->outmax[1], p->outmax[2]);
3893
3894 /* Do a specific test for out of Lab encoding range RGB primaries */
3895 /* (A more general check seems to get false positives - why ??) */
3896 if (h->pcs == icSigLabData
3897 && ( h->deviceClass == icSigDisplayClass
3898 || h->deviceClass == icSigOutputClass)
3899 && h->colorSpace == icSigRgbData) {
3900 double dev[3] = { 0.0, 0.0, 0.0 };
3901 double pcs[3];
3902 double clip = 0;
3903
3904 for (i = 0; i < 3; i++) {
3905 dev[i] = 1.0;
3906
3907 if (p->clut(p, pcs, dev) > 1)
3908 error ("%d, %s",p->pp->errc,p->pp->err);
3909
3910 /* Convert from efective pcs to natural pcs */
3911 if (p->pcs != icSigLabData)
3912 icmXYZ2Lab(&icmD50, pcs, pcs);
3913
3914 if (pcs[1] < -128.0 || pcs[1] > 128.0
3915 || pcs[2] < -128.0 || pcs[2] > 128.0) {
3916 warning("\n *** %s primary value can't be encoded in L*a*b* PCS (%f %f %f)",
3917 i == 0 ? "Red" : i == 1 ? "Green" : "Blue", pcs[0],pcs[1],pcs[2]);
3918 clip = 1;
3919 }
3920 dev[i] = 0.0;
3921 }
3922 if (clip)
3923 a1logw(g_log," *** Try switching to XYZ PCS ***\n");
3924 }
3925
3926 /* Use our rspl's to set the icc Lut AtoB table values. */
3927 /* Use helper function to do the hard work. */
3928 if (p->lut->set_tables(p->lut, ICM_CLUT_SET_EXACT, (void *)p,
3929 h->colorSpace, /* Input color space */
3930 h->pcs, /* Output color space */
3931 set_input, /* Input transfer function, Dev->Dev' */
3932 NULL, NULL, /* Use default Maximum range of Dev' values */
3933 set_clut, /* Dev' -> PCS' transfer function */
3934 NULL, NULL, /* Use default Maximum range of PCS' values */
3935 set_output, /* Linear output transform PCS'->PCS */
3936 NULL, NULL /* No APXLS */
3937 ) != 0)
3938 error("Setting 16 bit %s->%s Lut failed: %d, %s",
3939 icm2str(icmColorSpaceSignature, h->colorSpace),
3940 icm2str(icmColorSpaceSignature, h->pcs),
3941 p->pp->pp->errc,p->pp->pp->err);
3942
3943 #ifdef WARN_CLUT_CLIPPING
3944 if (p->pp->pp->warnc) {
3945 warning("Values clipped in setting A2B LUT!");
3946 if (p->pp->pp->warnc == 2 && h->pcs == icSigLabData) {
3947 a1logw(g_log,"*** Try switching to XYZ PCS ***\n");
3948 }
3949 }
3950 #endif /* WARN_CLUT_CLIPPING */
3951
3952 /* Init a CAM model in case it will be used (ie. in profile with camclip flag) */
3953 if (vc != NULL) /* One has been provided */
3954 p->vc = *vc; /* Copy the structure */
3955 else
3956 xicc_enum_viewcond(xicp, &p->vc, -1, NULL, 0, NULL); /* Use a default */
3957 p->cam = new_icxcam(cam_default);
3958 p->cam->set_view(p->cam, p->vc.Ev, p->vc.Wxyz, p->vc.La, p->vc.Yb, p->vc.Lv,
3959 p->vc.Yf, p->vc.Yg, p->vc.Gxyz, XICC_USE_HK, p->vc.hkscale);
3960
3961 if (flags & ICX_VERBOSE)
3962 printf("Done A to B table creation\n");
3963
3964 return (icxLuBase *)p;
3965 }
3966
3967 /* ========================================================== */
3968 /* Gamut boundary code. */
3969 /* ========================================================== */
3970
3971
3972 /* Context for creating gamut boundary points from, xicc */
3973 typedef struct {
3974 gamut *g; /* Gamut being created */
3975 icxLuLut *x; /* xLut we are working from */
3976 icxLuBase *flu; /* Forward xlookup */
3977 double in[MAX_CHAN]; /* Device input values */
3978 } lutgamctx;
3979
3980 /* Function to hand to zbrent to find a clut input' value at the ink limit */
3981 /* Returns value < 0.0 when within gamut, > 0.0 when out of gamut */
icxLimitFind(void * fdata,double tp)3982 static double icxLimitFind(void *fdata, double tp) {
3983 int i;
3984 double in[MAX_CHAN];
3985 lutgamctx *p = (lutgamctx *)fdata;
3986 double tt;
3987
3988 for (i = 0; i < p->x->inputChan; i++) {
3989 in[i] = tp * p->in[i]; /* Scale given input value */
3990 }
3991
3992 tt = icxLimitD(p->x, in);
3993
3994 return tt;
3995 }
3996
3997 /* Function to pass to rspl to create gamut boundary from */
3998 /* forward xLut transform grid points. */
3999 static void
lutfwdgam_func(void * pp,double * out,double * in)4000 lutfwdgam_func(
4001 void *pp, /* lutgamctx structure */
4002 double *out, /* output' value at clut grid point (ie. PCS' value) */
4003 double *in /* input' value at clut grid point (ie. device' value) */
4004 ) {
4005 lutgamctx *p = (lutgamctx *)pp;
4006 double pcso[3]; /* PCS output value */
4007
4008 /* Figure if we are over the ink limit. */
4009 if ( (p->x->ink.tlimit >= 0.0 || p->x->ink.klimit >= 0.0)
4010 && icxLimitD(p->x, in) > 0.0) {
4011 int i;
4012 double sf;
4013
4014 /* We are, so use the bracket search to discover a scale */
4015 /* for the clut input' value that will put us on the ink limit. */
4016
4017 for (i = 0; i < p->x->inputChan; i++)
4018 p->in[i] = in[i];
4019
4020 if (zbrent(&sf, 0.0, 1.0, 1e-4, icxLimitFind, pp) != 0) {
4021 return; /* Give up */
4022 }
4023
4024 /* Compute ink limit value */
4025 for (i = 0; i < p->x->inputChan; i++)
4026 p->in[i] = sf * in[i];
4027
4028 /* Compute the clut output for this clut input */
4029 p->x->clut(p->x, pcso, p->in);
4030 p->x->output(p->x, pcso, pcso);
4031 p->x->out_abs(p->x, pcso, pcso);
4032 } else { /* No ink limiting */
4033 /* Convert the clut PCS' values to PCS output values */
4034 p->x->output(p->x, pcso, out);
4035 p->x->out_abs(p->x, pcso, pcso);
4036 }
4037
4038 /* Expand the gamut surface with this point */
4039 p->g->expand(p->g, pcso);
4040
4041 /* Leave out[] unchanged */
4042 }
4043
4044
4045 /* Function to pass to rspl to create gamut boundary from */
4046 /* backwards Lut transform. This is called for every node in the */
4047 /* B2A grid. */
4048 static void
lutbwdgam_func(void * pp,double * out,double * in)4049 lutbwdgam_func(
4050 void *pp, /* lutgamctx structure */
4051 double *out, /* output value */
4052 double *in /* input value */
4053 ) {
4054 lutgamctx *p = (lutgamctx *)pp;
4055 double devo[MAX_CHAN]; /* Device output value */
4056 double pcso[3]; /* PCS output value */
4057
4058 /* Convert the clut values to device output values */
4059 p->x->output(p->x, devo, out); /* (Device never uses out_abs()) */
4060
4061 /* Convert from device values to PCS values */
4062 p->flu->lookup(p->flu, pcso, devo);
4063
4064 /* Expand the gamut surface with this point */
4065 p->g->expand(p->g, pcso);
4066
4067 /* Leave out[] unchanged */
4068 }
4069
4070 /* Given an xicc lookup object, return a gamut object. */
4071 /* Note that the PCS must be Lab or Jab */
4072 /* An icxLuLut type must be icmFwd or icmBwd, */
4073 /* and for icmFwd, the ink limit (if supplied) will be applied. */
4074 /* Return NULL on error, check errc+err for reason */
icxLuLutGamut(icxLuBase * plu,double detail)4075 static gamut *icxLuLutGamut(
4076 icxLuBase *plu, /* this */
4077 double detail /* gamut detail level, 0.0 = def */
4078 ) {
4079 xicc *p = plu->pp; /* parent xicc */
4080 icxLuLut *luluto = (icxLuLut *)plu; /* Lookup xLut type object */
4081 icColorSpaceSignature ins, pcs, outs;
4082 icmLookupFunc func;
4083 icRenderingIntent intent;
4084 double white[3], black[3], kblack[3];
4085 int inn, outn;
4086 gamut *gam;
4087
4088 /* get some details */
4089 plu->spaces(plu, &ins, &inn, &outs, &outn, NULL, &intent, &func, &pcs);
4090
4091 if (func != icmFwd && func != icmBwd) {
4092 p->errc = 1;
4093 sprintf(p->err,"Creating Gamut surface for anything other than Device <-> PCS is not supported.");
4094 return NULL;
4095 }
4096
4097 if (pcs != icSigLabData && pcs != icxSigJabData) {
4098 p->errc = 1;
4099 sprintf(p->err,"Creating Gamut surface PCS of other than Lab or Jab is not supported.");
4100 return NULL;
4101 }
4102
4103 if (func == icmFwd) {
4104 lutgamctx cx;
4105
4106 cx.g = gam = new_gamut(detail, pcs == icxSigJabData, 0);
4107 cx.x = luluto;
4108
4109 /* Scan through grid. */
4110 /* (Note this can give problems for a strange input space - ie. Lab */
4111 /* and a low grid resolution - ie. 2) */
4112 luluto->clutTable->scan_rspl(
4113 luluto->clutTable, /* this */
4114 RSPL_NOFLAGS, /* Combination of flags */
4115 (void *)&cx, /* Opaque function context */
4116 lutfwdgam_func /* Function to set from */
4117 );
4118
4119 /* Make sure the white and point goes in too, if it isn't in the grid */
4120 plu->efv_wh_bk_points(plu, white, NULL, NULL);
4121 gam->expand(gam, white);
4122
4123 if (detail == 0.0)
4124 detail = 10.0;
4125
4126 /* If the gamut is more than cursary, add some more detail surface points */
4127 if (detail < 20.0 || luluto->clutTable->g.mres < 4) {
4128 int res;
4129 DCOUNT(co, MAX_CHAN, inn, 0, 0, 2);
4130
4131 res = (int)(500.0/detail); /* Establish an appropriate sampling density */
4132 if (res < 10)
4133 res = 10;
4134
4135 /* Itterate over all the faces in the device space */
4136 DC_INIT(co);
4137 while(!DC_DONE(co)) { /* Count through the corners of hyper cube */
4138 int e, m1, m2;
4139 double in[MAX_CHAN];
4140 double out[3];
4141
4142 for (e = 0; e < inn; e++) { /* Base value */
4143 in[e] = (double)co[e]; /* Base value */
4144 in[e] = in[e] * (luluto->inmax[e] - luluto->inmin[e])
4145 + luluto->inmin[e];
4146 }
4147
4148 /* Figure if we are over the ink limit. */
4149 if ((luluto->ink.tlimit >= 0.0 || luluto->ink.klimit >= 0.0)
4150 && icxLimit(luluto, in) > 0.0) {
4151 DC_INC(co);
4152 continue; /* Skip points over limit */
4153 }
4154
4155 /* Scan only device surface */
4156 for (m1 = 0; m1 < (inn-1); m1++) { /* Choose first coord to scan */
4157 if (co[m1] != 0)
4158 continue; /* Not at lower corner */
4159 for (m2 = m1 + 1; m2 < inn; m2++) { /* Choose second coord to scan */
4160 int x, y;
4161
4162 if (co[m2] != 0)
4163 continue; /* Not at lower corner */
4164
4165 for (x = 0; x < res; x++) { /* step over surface */
4166 in[m1] = x/(res - 1.0);
4167 in[m1] = in[m1] * (luluto->inmax[m1] - luluto->inmin[m1])
4168 + luluto->inmin[m1];
4169 for (y = 0; y < res; y++) {
4170 in[m2] = y/(res - 1.0);
4171 in[m2] = in[m2] * (luluto->inmax[m2] - luluto->inmin[m2])
4172 + luluto->inmin[m2];
4173
4174 /* Figure if we are over the ink limit. */
4175 if ( (luluto->ink.tlimit >= 0.0 || luluto->ink.klimit >= 0.0)
4176 && icxLimit(luluto, in) > 0.0) {
4177 continue; /* Skip points over limit */
4178 }
4179
4180 luluto->lookup((icxLuBase *)luluto, out, in);
4181 gam->expand(gam, out);
4182 }
4183 }
4184 }
4185 }
4186 /* Increment index within block */
4187 DC_INC(co);
4188 }
4189 }
4190
4191 /* Now set the cusp points by itterating through colorant 0 & 100% combinations */
4192 /* If we know what sort of space it is: */
4193 if (ins == icSigRgbData || ins == icSigCmyData || ins == icSigCmykData) {
4194 DCOUNT(co, 3, 3, 0, 0, 2);
4195
4196 gam->setcusps(gam, 0, NULL);
4197 DC_INIT(co);
4198 while(!DC_DONE(co)) {
4199 int e;
4200 double in[MAX_CHAN];
4201 double out[3];
4202
4203 if (!(co[0] == 0 && co[1] == 0 && co[2] == 0)
4204 && !(co[0] == 1 && co[1] == 1 && co[2] == 1)) { /* Skip white and black */
4205 for (e = 0; e < 3; e++)
4206 in[e] = (double)co[e];
4207 in[e] = 0; /* K */
4208
4209 /* Always use the device->PCS conversion */
4210 if (luluto->lookup((icxLuBase *)luluto, out, in) > 1)
4211 error ("%d, %s",p->errc,p->err);
4212 gam->setcusps(gam, 3, out);
4213 }
4214
4215 DC_INC(co);
4216 }
4217 gam->setcusps(gam, 2, NULL);
4218
4219 /* Do all ink combinations and hope we can sort it out */
4220 /* (This may not be smart, since bodgy cusp values will cause gamut mapping to fail...) */
4221 } else if (ins != icSigXYZData
4222 && ins != icSigLabData
4223 && ins != icSigLuvData
4224 && ins != icSigYxyData) {
4225 DCOUNT(co, MAX_CHAN, inn, 0, 0, 2);
4226
4227 gam->setcusps(gam, 0, NULL);
4228 DC_INIT(co);
4229 while(!DC_DONE(co)) {
4230 int e;
4231 double in[MAX_CHAN];
4232 double out[3];
4233
4234 for (e = 0; e < inn; e++) {
4235 in[e] = (double)co[e];
4236 in[e] = in[e] * (luluto->inmax[e] - luluto->inmin[e])
4237 + luluto->inmin[e];
4238 }
4239
4240 /* Figure if we are over the ink limit. */
4241 if ((luluto->ink.tlimit >= 0.0 || luluto->ink.klimit >= 0.0)
4242 && icxLimit(luluto, in) > 0.0) {
4243 DC_INC(co);
4244 continue; /* Skip points over limit */
4245 }
4246
4247 luluto->lookup((icxLuBase *)luluto, out, in);
4248 gam->setcusps(gam, 1, out);
4249
4250 DC_INC(co);
4251 }
4252 gam->setcusps(gam, 2, NULL);
4253 }
4254
4255 } else { /* Must be icmBwd */
4256 lutgamctx cx;
4257
4258 /* Get an appropriate device to PCS conversion for the fwd conversion */
4259 /* we use after bwd conversion in lutbwdgam_func() */
4260 switch (intent) {
4261 /* If it is relative */
4262 case icmDefaultIntent: /* Shouldn't happen */
4263 case icPerceptual:
4264 case icRelativeColorimetric:
4265 case icSaturation:
4266 intent = icRelativeColorimetric; /* Choose relative */
4267 break;
4268 /* If it is absolute */
4269 case icAbsoluteColorimetric:
4270 case icxAppearance:
4271 case icxAbsAppearance:
4272 break; /* Leave unchanged */
4273 default:
4274 break;
4275 }
4276 if ((cx.flu = p->get_luobj(p, ICX_CLIP_NEAREST, icmFwd, intent, pcs, icmLuOrdNorm,
4277 &plu->vc, NULL)) == NULL) {
4278 return NULL; /* oops */
4279 }
4280
4281 cx.g = gam = new_gamut(detail, pcs == icxSigJabData, 0);
4282
4283 cx.x = luluto;
4284
4285 luluto->clutTable->scan_rspl(
4286 luluto->clutTable, /* this */
4287 RSPL_NOFLAGS, /* Combination of flags */
4288 (void *)&cx, /* Opaque function context */
4289 lutbwdgam_func /* Function to set from */
4290 );
4291
4292 /* Now set the cusp points by using the fwd conversion and */
4293 /* itterating through colorant 0 & 100% combinations. */
4294 /* If we know what sort of space it is: */
4295 if (outs == icSigRgbData || outs == icSigCmyData || outs == icSigCmykData) {
4296 DCOUNT(co, 3, 3, 0, 0, 2);
4297
4298 gam->setcusps(gam, 0, NULL);
4299 DC_INIT(co);
4300 while(!DC_DONE(co)) {
4301 int e;
4302 double in[MAX_CHAN];
4303 double out[3];
4304
4305 if (!(co[0] == 0 && co[1] == 0 && co[2] == 0)
4306 && !(co[0] == 1 && co[1] == 1 && co[2] == 1)) { /* Skip white and black */
4307 for (e = 0; e < 3; e++)
4308 in[e] = (double)co[e];
4309 in[e] = 0; /* K */
4310
4311 /* Always use the device->PCS conversion */
4312 if (cx.flu->lookup((icxLuBase *)cx.flu, out, in) > 1)
4313 error ("%d, %s",p->errc,p->err);
4314 gam->setcusps(gam, 3, out);
4315 }
4316
4317 DC_INC(co);
4318 }
4319 gam->setcusps(gam, 2, NULL);
4320
4321 /* Do all ink combinations and hope we can sort it out */
4322 /* (This may not be smart, since bodgy cusp values will cause gamut mapping to fail...) */
4323 } else if (ins != icSigXYZData
4324 && ins != icSigLabData
4325 && ins != icSigLuvData
4326 && ins != icSigYxyData) {
4327 DCOUNT(co, MAX_CHAN, inn, 0, 0, 2);
4328
4329 gam->setcusps(gam, 0, NULL);
4330 DC_INIT(co);
4331 while(!DC_DONE(co)) {
4332 int e;
4333 double in[MAX_CHAN];
4334 double out[3];
4335
4336 for (e = 0; e < inn; e++) {
4337 in[e] = (double)co[e];
4338 in[e] = in[e] * (luluto->inmax[e] - luluto->inmin[e])
4339 + luluto->inmin[e];
4340 }
4341
4342 /* Figure if we are over the ink limit. */
4343 if ((luluto->ink.tlimit >= 0.0 || luluto->ink.klimit >= 0.0)
4344 && icxLimit(luluto, in) > 0.0) {
4345 DC_INC(co);
4346 continue; /* Skip points over limit */
4347 }
4348
4349 cx.flu->lookup((icxLuBase *)cx.flu, out, in);
4350 gam->setcusps(gam, 1, out);
4351
4352 DC_INC(co);
4353 }
4354 gam->setcusps(gam, 2, NULL);
4355 }
4356 cx.flu->del(cx.flu); /* Done with the fwd conversion */
4357 }
4358
4359 /* Set the white and black points */
4360 plu->efv_wh_bk_points(plu, white, black, kblack);
4361 gam->setwb(gam, white, black, kblack);
4362
4363 //printf("~1 icxLuLutGamut: set black %f %f %f and kblack %f %f %f\n", black[0], black[1], black[2], kblack[0], kblack[1], kblack[2]);
4364
4365 #ifdef NEVER /* Not sure if this is a good idea ?? */
4366 /* Since we know this is a profile, force the space and gamut points to be the same */
4367 gam->getwb(gam, NULL, NULL, white, black, kblack); /* Get the actual gamut white and black points */
4368 gam->setwb(gam, white, black, kblack); /* Put it back as colorspace one */
4369 #endif
4370
4371 return gam;
4372 }
4373
4374 /* ========================================================== */
4375 /* Cusp Map finding code. */
4376 /* ========================================================== */
4377
4378 /* Context for creating Cusp Map points from, xicc */
4379 typedef struct {
4380 icxCuspMap *cm; /* Cusp Map being created */
4381 icxLuLut *x; /* xLut we are working from */
4382 icxLuBase *flu; /* Forward xlookup */
4383 double in[MAX_CHAN]; /* Device input values */
4384 } lutcuspmapctx;
4385
4386 /* Function to pass to rspl to create cusp map from */
4387 /* forward xLut transform grid points. */
4388 static void
lutfwdcuspmap_func(void * pp,double * out,double * in)4389 lutfwdcuspmap_func(
4390 void *pp, /* lutcuspmapctx structure */
4391 double *out, /* output' value at clut grid point (ie. PCS' value) */
4392 double *in /* input' value at clut grid point (ie. device' value) */
4393 ) {
4394 lutcuspmapctx *p = (lutcuspmapctx *)pp;
4395 double pcso[3]; /* PCS output value */
4396
4397 /* Figure if we are over the ink limit. */
4398 if ( (p->x->ink.tlimit >= 0.0 || p->x->ink.klimit >= 0.0)
4399 && icxLimitD(p->x, in) > 0.0) {
4400 int i;
4401 double sf;
4402
4403 /* We are, so use the bracket search to discover a scale */
4404 /* for the clut input' value that will put us on the ink limit. */
4405
4406 for (i = 0; i < p->x->inputChan; i++)
4407 p->in[i] = in[i];
4408
4409 if (zbrent(&sf, 0.0, 1.0, 1e-4, icxLimitFind, pp) != 0) {
4410 return; /* Give up */
4411 }
4412
4413 /* Compute ink limit value */
4414 for (i = 0; i < p->x->inputChan; i++)
4415 p->in[i] = sf * in[i];
4416
4417 /* Compute the clut output for this clut input */
4418 p->x->clut(p->x, pcso, p->in);
4419 p->x->output(p->x, pcso, pcso);
4420 p->x->out_abs(p->x, pcso, pcso);
4421 } else { /* No ink limiting */
4422 /* Convert the clut PCS' values to PCS output values */
4423 p->x->output(p->x, pcso, out);
4424 p->x->out_abs(p->x, pcso, pcso);
4425 }
4426
4427 /* Expand the usp map surface with this point */
4428 p->cm->expand(p->cm, pcso);
4429
4430 /* Leave out[] unchanged */
4431 }
4432
4433 /* Expand cusp map with given point */
4434 /* (We use a segmentned maxima approach) */
cuspmap_expand(icxCuspMap * p,double lab[3])4435 static void cuspmap_expand(icxCuspMap *p, double lab[3]) {
4436 double h, C;
4437 int ix;
4438
4439 /* Hue angle 0.0 .. 1.0 */
4440 h = (0.5/3.14159265359) * atan2(lab[2], lab[1]);
4441 h = (h < 0.0) ? h + 1.0 : h;
4442
4443 /* Slot index 0..res-1 */
4444 ix = (int)floor(h * p->res + 0.5);
4445 if (ix >= p->res)
4446 ix -= p->res;
4447
4448 /* Chromanance */
4449 C = sqrt(lab[1] * lab[1] + lab[2] * lab[2]);
4450
4451 /* New point at this angle with largest chromanance */
4452 if (C > p->C[ix]) {
4453 p->C[ix] = C;
4454 p->L[ix] = lab[0];
4455 }
4456
4457 /* Tracl min * max L */
4458 if (lab[0] > p->Lmax[0])
4459 icmCpy3(p->Lmax, lab);
4460 if (lab[0] < p->Lmin[0])
4461 icmCpy3(p->Lmin, lab);
4462 }
4463
4464 /* Interpolate over any gaps in map */
cuspmap_complete(icxCuspMap * p)4465 static void cuspmap_complete(icxCuspMap *p) {
4466 int i, j, k;
4467
4468 // printf("cuspmap list before fixups:\n");
4469 // for (i = 0; i < p->res; i++) {
4470 // printf(" %d: C = %f, L = %f\n",i, p->C[i], p->L[i]);
4471 // }
4472
4473 /* First check if there are any entries at all */
4474 for (i = 0; i < p->res; i++) {
4475 if (p->C[i] >= 0.0)
4476 break;
4477 }
4478
4479 if (i >= p->res) /* Nothing there - give up */
4480 return;
4481
4482 /* See if there are any gaps */
4483 j = -1;
4484 for (i = 0; i < p->res; i++) {
4485 //printf("~1 checking %d\n",i);
4486 if (p->C[i] >= 0.0) {
4487 j = i; /* Last valid slot */
4488 //printf("~1 last valid %d\n",j);
4489
4490 /* Got a gap */
4491 } else {
4492 int ii = i;
4493 //printf("~1 found gap at %d\n",i);
4494 if (j < 0) { /* No previous */
4495 //printf("~1 no previous\n");
4496 for (j = p->res-1; j >= 0; j--) {
4497 if (p->C[j] >= 0.0)
4498 break;
4499 }
4500 }
4501 //printf("~1 got previous %d\n",j);
4502 /* Find next, even if it's behind us */
4503 for (k = i+1; k != i; k = (k+1) % p->res) {
4504 if (p->C[k] >= 0.0)
4505 break;
4506 }
4507 //printf("~1 got next %d\n",k);
4508
4509 /* Interpolate between them */
4510 for (; i != k; i = (i+1) % p->res) {
4511 double prop;
4512 int bb, tt;
4513 bb = k > i ? k - i : k + p->res - i;
4514 tt = k > j ? k - j : k + p->res - j;
4515 prop = (double)bb/(double)tt;
4516 p->C[i] = prop * p->C[j] + (1.0 - prop) * p->C[k];
4517 p->L[i] = prop * p->L[j] + (1.0 - prop) * p->L[k];
4518 //printf("~1 interpolating %d with weight %f from %d and %d\n",i,prop,j,k);
4519 }
4520 if (k > ii) {
4521 i = k; /* Continue from next valid */
4522 //printf("~1 continuing from %d\n",i);
4523 } else {
4524 //printf("~1 we've looped back\n");
4525 break; /* We've looped back */
4526 }
4527 }
4528 }
4529
4530 // printf("cuspmap list after fixups:\n");
4531 // for (i = 0; i < p->res; i++) {
4532 // printf(" %d: C = %f, L = %f\n",i, p->C[i], p->L[i]);
4533 // }
4534 }
4535
4536 /* Return the corresponding cusp location, given the source point */
cuspmap_getCusp(icxCuspMap * p,double cuspLCh[3],double srcLab[3])4537 static void cuspmap_getCusp(icxCuspMap *p,double cuspLCh[3], double srcLab[3]) {
4538 double h, C;
4539 int ix, ix0, ix1;
4540
4541 /* Hue angle 0.0 .. 1.0 */
4542 h = (0.5/3.14159265359) * atan2(srcLab[2], srcLab[1]);
4543 h = (h < 0.0) ? h + 1.0 : h;
4544
4545 //printf("~1 getCusp in %f %f %f, h = %f\n", srcLab[0], srcLab[1], srcLab[2],h);
4546
4547 /* Slot index 0..res-1 */
4548 ix = (int)floor(h * p->res + 0.5);
4549 if (ix >= p->res)
4550 ix -= p->res;
4551
4552 /* Indexes each side */
4553 ix0 = ix > 0 ? ix-1 : p->res-1;
4554 ix1 = ix < (p->res-1) ? ix+1 : 0;
4555
4556 cuspLCh[0] = p->L[ix];
4557 cuspLCh[1] = p->C[ix];
4558 if (cuspLCh[1] > p->C[ix0]) /* Be conservative with C value */
4559 cuspLCh[1] = p->C[ix0];
4560 if (cuspLCh[1] > p->C[ix1])
4561 cuspLCh[1] = p->C[ix1];
4562 cuspLCh[2] = 360.0 * h;
4563
4564 //printf("~1 index %d, L %f C %f h %f\n", ix, cuspLCh[0], cuspLCh[1], cuspLCh[2]);
4565 }
4566
4567 /* We're done with CuspMap */
cuspmap_del(icxCuspMap * p)4568 static void cuspmap_del(icxCuspMap *p) {
4569 if (p != NULL) {
4570 if (p->L != NULL)
4571 free(p->L);
4572 if (p->C != NULL)
4573 free(p->C);
4574 free(p);
4575 }
4576 }
4577
4578 /* Given an xicc lookup object, return an icxCuspMap object. */
4579 /* Note that the PCS must be Lab or Jab. */
4580 /* An icxLuLut type must be icmFwd, and the ink limit (if supplied) */
4581 /* will be applied. */
4582 /* Return NULL on error, check errc+err for reason */
icxLuLutCuspMap(icxLuBase * plu,int res)4583 static icxCuspMap *icxLuLutCuspMap(
4584 icxLuBase *plu, /* this */
4585 int res /* Hue resolution */
4586 ) {
4587 xicc *p = plu->pp; /* parent xicc */
4588 icxLuLut *luluto = (icxLuLut *)plu; /* Lookup xLut type object */
4589 icColorSpaceSignature ins, pcs, outs;
4590 icmLookupFunc func;
4591 icRenderingIntent intent;
4592 int inn, outn;
4593 lutcuspmapctx cx;
4594 icxCuspMap *cm;
4595 int i;
4596
4597 /* get some details */
4598 plu->spaces(plu, &ins, &inn, &outs, &outn, NULL, &intent, &func, &pcs);
4599
4600 if (func != icmFwd) {
4601 p->errc = 1;
4602 sprintf(p->err,"Creating CuspMap for anything other than Device -> PCS is not supported.");
4603 return NULL;
4604 }
4605
4606 if (pcs != icSigLabData && pcs != icxSigJabData) {
4607 p->errc = 1;
4608 sprintf(p->err,"Creating CuspMap PCS of other than Lab or Jab is not supported.");
4609 return NULL;
4610 }
4611
4612 cx.cm = cm = (icxCuspMap *) calloc(1, sizeof(icxCuspMap));
4613 cx.x = luluto;
4614
4615 if (cx.cm == NULL) {
4616 p->errc = 2;
4617 sprintf(p->err,"Malloc of icxCuspMap failed");
4618 return NULL;
4619 }
4620
4621 if ((cm->L = (double *)malloc(sizeof(double) * res)) == NULL) {
4622 free(cm);
4623 p->errc = 2;
4624 sprintf(p->err,"Malloc of icxCuspMap failed");
4625 return NULL;
4626 }
4627
4628 if ((cm->C = (double *)malloc(sizeof(double) * res)) == NULL) {
4629 free(cm->L);
4630 free(cm);
4631 p->errc = 2;
4632 sprintf(p->err,"Malloc of icxCuspMap failed");
4633 return NULL;
4634 }
4635
4636 cm->res = res;
4637 cm->expand = cuspmap_expand;
4638 cm->getCusp = cuspmap_getCusp;
4639 cm->del = cuspmap_del;
4640
4641 for (i = 0; i < res; i++) {
4642 cm->L[i] = -1.0;
4643 cm->C[i] = -1.0;
4644 }
4645 cm->Lmax[0] = -1.0;
4646 cm->Lmin[0] = 101.0;
4647
4648 /* Scan through grid, expanding the CuspMap. */
4649 luluto->clutTable->scan_rspl(
4650 luluto->clutTable, /* this */
4651 RSPL_NOFLAGS, /* Combination of flags */
4652 (void *)&cx, /* Opaque function context */
4653 lutfwdcuspmap_func /* Function to set from */
4654 );
4655
4656 cm->Lmax[0] -= 0.1; /* Make sure tips intersect */
4657 cm->Lmin[0] += 0.1;
4658
4659 for (i = 0; i < res; i++) {
4660 if (cm->L[i] > cm->Lmax[0])
4661 cm->L[i] = cm->Lmax[0];
4662 if (cm->L[i] < cm->Lmin[0])
4663 cm->L[i] = cm->Lmin[0];
4664 }
4665
4666 /* Fill in any gaps */
4667 cuspmap_complete(cm);
4668
4669 return cm;
4670 }
4671
4672 /* ----------------------------------------------- */
4673 #ifdef DEBUG
4674 #undef DEBUG /* Limit extent to this file */
4675 #endif
4676
4677
4678
4679
4680
4681
4682