1
2
3 /*
4 * International Color Consortium color transform expanded support
5 *
6 * Author: Graeme W. Gill
7 * Date: 2/7/00
8 * Version: 1.00
9 *
10 * Copyright 2000 - 2007 Graeme W. Gill
11 * All rights reserved.
12 * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
13 * see the License.txt file for licencing details.
14 *
15 * Based on the xlut.c code.
16 */
17
18 /*
19 * TTBD:
20 *
21 * Need to use this for B2A tables rather than inverting
22 * A2B curves. Need to add grid sizing to cover just gamut range
23 * (including level axis gamut, but watch out for devices that
24 * have values below the black point or above the white point),
25 * 3x3 matrix optimization, and white point to grid node mapping for B2A.
26 * Currently code assumes output is always PCS ?? - would need to fix for opposite.
27 *
28 * Currently the Lab A2B output tables are adjusted for ab symetry
29 * to make the B2A white point land on a grid point, given that
30 * the icc code forces a symetric ab range. This only works
31 * because the B2A is using the same per channel curves.
32 * Done properly, it should be possible to know where grid
33 * points land within the range, and to modify the B2A input curves
34 * to make the white point land on a grid point.
35 * (For B2A the pseudo least squares adjustment needs to be turned
36 * off for that grid point too.)
37 *
38 * Note that one quandry is that the curve fitting doesn't
39 * fit well when the input data has an offset and/or plateaus.
40 */
41
42 /*
43 * This module provides curve and matrix fitting functionality used to
44 * create per channel input and/or output curves for clut profiles,
45 * and optionally creates the rspl to fit the data as well.
46 *
47 * The approach used is to initialy creating input and output shaper
48 * curves that minimize the overall delta E of the test point set,
49 * when a linear matrix is substituted for a clut is used.
50 * (This is the same as pre-V0.70 approach. )
51 *
52 * The residual error between the test point set and the shaper/matrix/shaper
53 * model is then computed, and mapped against each input channel, and
54 * a positioning mapping curve created that aims to map rspl grid
55 * locations more densly where residual errors are high, and more
56 * sparsely where they are low.
57 * The input shaper and positioning curves are then combined together,
58 * so that the positioning curve determines which input values map
59 * to rspl grid points, and the shaper curve determines the mapping
60 * between the grid points. The width of the grid cells is computed
61 * in the shaper mapped input space, and then fed to the rspl fitting
62 * code, so that its smoothness evaluation function can take account
63 * of the non-uniform spacing of the grid points.
64 */
65
66 #include <sys/types.h>
67 #include <string.h>
68 #include <ctype.h>
69 #ifdef __sun
70 #include <unistd.h>
71 #endif
72 #include "copyright.h"
73 #include "aconfig.h"
74 #include "numlib.h"
75 #include "icc.h"
76 #include "rspl.h"
77 #include "xicc.h"
78 #include "plot.h"
79 #include "xfit.h"
80 #include "sort.h"
81
82 #undef USE_XYZ_Y2LCURVE /* [Und] Use underlying L* curve for XYZ encoding */
83 /* This seems to work badly, even with high smoothness. Why ? */
84 /* It does speed up 1D lut creation though. */
85
86 #undef DEBUG /* Debug information */
87 #undef DEBUG_PROGRESS /* Show powell progress */
88 #undef DEBUG_PLOT /* Plot in & out curves */
89 #undef SPECIAL_FORCE /* Check rspl nodes against linear XYZ model */
90 #undef SPECIAL_FORCE_GAMMA /* Force correct gamma shaper curves */
91 #undef SPECIAL_TEST_GAMMA /* Use gamma on reference model */
92 #undef SPECIAL_TEST_LAB
93
94 #define RSPLFLAGS (0 /* | RSPL_2PASSSMTH */ /* | RSPL_EXTRAFIT2 */)
95
96 #undef EXTEND_GRID /* [Undef] Use extended RSPL grid around interpolation */
97 #define EXTEND_GRID_BYN 2 /* Rows to extend grid. */
98
99 #undef NODDV /* Use slow non d/dv powell else use conjgrad */
100 #define CURVEPOW 1.0 /* Power to raise deltaE squared to in setting in/out curves */
101 /* This provides a means of punishing high maximum errors. */
102
103 #define POWTOL1 1e-3 /* Shaper Powell optimiser tollerance for first passes */
104 #define MAXITS1 1000 /* Shaper number of itterations for first passes */
105 #define POWTOL 1e-5 /* Shaper Powell optimiser tollerance in delta E squared ^ CURVEPOW */
106 #define MAXITS 4000 /* Shaper number of itterations before giving up */
107 #define PDDEL 1e-6 /* Fake partial derivative del */
108
109 /* Weights for shaper in/out curve parameters, to minimise unconstrained "wiggles" */
110 #define SHAPE_WEIGHT 1.0 /* Overal shaper weight contribution - err on side of smoothness */
111 #define SHAPE_HW01 0.002 /* 0 & 1 harmonic weights */
112 #define SHAPE_HBREAK 4 /* Harmonic that has HWBR */
113 #define SHAPE_HWBR 20.0 /* Base weight of harmonics HBREAK up */
114 #define SHAPE_HWINC 60.0 /* Increase in weight for each harmonic above HWBR */
115
116 /* Weights for the positioning curve parameters */
117 #define PSHAPE_MINE 0.02 /* Minum background residual error level */
118 #define PSHAPE_DIST 1.0 /* Agressivness of grid distribution */
119
120
121 /* - - - - - - - - - - - - - - - - - */
122
123 /* Extra non-linearity used as base for XYZ output curves. */
124 /* This makes the XYZ grid values more perceptual, and asks less */
125 /* of the automatically created output curve shape. */
126 /* (We assume XYZ is in 0..1 scale */
127
128 /* Transfer function with offset and scale + Y2L curve */
icxSTransFuncY2L(double * v,int luord,double vv,double min,double max)129 static double icxSTransFuncY2L(
130 double *v, /* Pointer to first parameter */
131 int luord, /* Number of parameters */
132 double vv, /* Source of value */
133 double min, /* Scale values */
134 double max
135 ) {
136 max -= min;
137
138 vv = (vv - min)/max;
139
140 #ifdef USE_XYZ_Y2LCURVE
141 if (vv > 0.008856451586)
142 vv = 1.16 * pow(vv,1.0/3.0) - 0.16;
143 else
144 vv = 9.032962896 * vv;
145 #endif
146
147 vv = icxTransFunc(v, luord, vv);
148
149 vv = (vv * max) + min;
150 return vv;
151 }
152
153 /* Inverse Transfer function with offset and scale + Y2L */
icxInvSTransFuncY2L(double * v,int luord,double vv,double min,double max)154 static double icxInvSTransFuncY2L(
155 double *v, /* Pointer to first parameter */
156 int luord, /* Number of parameters */
157 double vv, /* Source of value */
158 double min, /* Scale values */
159 double max
160 ) {
161 max -= min;
162
163 vv = (vv - min)/max;
164 vv = icxInvTransFunc(v, luord, vv);
165
166 #ifdef USE_XYZ_Y2LCURVE
167 if (vv > 0.08)
168 vv = pow((vv + 0.16)/1.16, 3.0);
169 else
170 vv = vv/9.032962896;
171 #endif
172
173 vv = (vv * max) + min;
174
175 return vv;
176 }
177
178 /* Transfer function with offset and scale, and */
179 /* partial derivative with respect to the */
180 /* parameters and the input value. */
icxdpdiSTransFuncY2L(double * v,double * dv,double * pdin,int luord,double vv,double min,double max)181 static double icxdpdiSTransFuncY2L(
182 double *v, /* Pointer to first parameter */
183 double *dv, /* Return derivative wrt each parameter */
184 double *pdin, /* Return derivative wrt source value */
185 int luord, /* Number of parameters */
186 double vv, /* Source of value */
187 double min, /* Scale values */
188 double max
189 ) {
190 int i;
191 double idv = 1.0;
192 max -= min;
193
194 #ifdef USE_XYZ_Y2LCURVE
195 if (vv > 0.008856451586) {
196 vv = 1.16 * pow(vv,1.0/3.0) - 0.16;
197 idv = 1.16 / (3.0 * pow(vv, 2.0/3.0));
198 } else {
199 vv = 9.032962896 * vv;
200 idv = 9.032962896;
201 }
202 #endif
203
204 vv = (vv - min)/max;
205
206 vv = icxdpdiTransFunc(v, dv, pdin, luord, vv);
207
208 *pdin *= idv; /* Account for input multiplier */
209
210 vv = (vv * max) + min;
211
212 for (i = 0; i < luord; i++) {
213 dv[i] *= max;
214 }
215 return vv;
216 }
217
218 /* - - - - - - - - - - - - - - - - - */
219
220 #ifdef DEBUG
221 static void dump_xfit(xfit *p);
222 #endif
223
224 #ifdef DEBUG_PLOT /* Not currently used in runtime code*/
225
226 /* Lookup a value though an input position curve */
xfit_poscurve(xfit * p,double in,int chan)227 static double xfit_poscurve(xfit *p, double in, int chan) {
228 double rv = in;
229 if (p->tcomb & oc_p)
230 rv = icxSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], rv,
231 p->in_min[chan], p->in_max[chan]);
232 return rv;
233 }
234
235 /* Lookup di values though the input position curves */
xfit_poscurves(xfit * p,double * out,double * in)236 static void xfit_poscurves(xfit *p, double *out, double *in) {
237 int e;
238
239 for (e = 0; e < p->di; e++) {
240 double val = in[e];
241 if (p->tcomb & oc_i)
242 val = icxSTransFunc(p->v + p->pos_offs[e], p->iluord[e], val,
243 p->in_min[e], p->in_max[e]);
244 out[e] = val;
245 }
246 }
247
248 /* Inverse Lookup a value though an input position curve */
xfit_invposcurve(xfit * p,double in,int chan)249 static double xfit_invposcurve(xfit *p, double in, int chan) {
250 double rv = in;
251 if (p->tcomb & oc_i)
252 rv = icxInvSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], rv,
253 p->in_min[chan], p->in_max[chan]);
254 return rv;
255 }
256
257 /* Inverse Lookup di values though input position curves */
xfit_invposcurves(xfit * p,double * out,double * in)258 static void xfit_invposcurves(xfit *p, double *out, double *in) {
259 int e;
260
261 for (e = 0; e < p->di; e++) {
262 double val = in[e];
263 if (p->tcomb & oc_i)
264 val = icxInvSTransFunc(p->v + p->pos_offs[e], p->iluord[e], val,
265 p->in_min[e], p->in_max[e]);
266 out[e] = val;
267 }
268 }
269 #endif /* DEBUG_PLOT */
270
271 /* - - - - - - - - - - - - - - - - - */
272 /* Lookup a value though input shape curve */
xfit_shpcurve(xfit * p,double in,int chan)273 static double xfit_shpcurve(xfit *p, double in, int chan) {
274 double rv = in;
275
276 if (p->tcomb & oc_i) {
277 #ifdef SPECIAL_FORCE_GAMMA
278 double gam;
279 if (chan == 0)
280 gam = 1.9;
281 else if (chan == 1)
282 gam = 2.0;
283 else if (chan == 2)
284 gam = 2.1;
285 else
286 gam = 1.0;
287 rv = pow(in, gam);
288 #else
289 rv = icxSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], rv,
290 p->in_min[chan], p->in_max[chan]);
291 #endif
292 }
293 return rv;
294 }
295
296 #ifdef NEVER /* Not currently used */
297 /* Lookup a value though shape curves */
xfit_shpcurves(xfit * p,double * out,double * in)298 static void xfit_shpcurves(xfit *p, double *out, double *in) {
299 int e;
300
301 for (e = 0; e < p->di; e++) {
302 double val = in[e];
303 if (p->tcomb & oc_i)
304 val = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], val,
305 p->in_min[e], p->in_max[e]);
306 out[e] = val;
307 }
308 }
309 #endif /* NEVER */
310
311 /* Inverse Lookup a value though a shape curve */
xfit_invshpcurve(xfit * p,double in,int chan)312 static double xfit_invshpcurve(xfit *p, double in, int chan) {
313 double rv = in;
314
315 if (p->tcomb & oc_i) {
316 #ifdef SPECIAL_FORCE_GAMMA
317 double gam;
318 if (chan == 0)
319 gam = 1.9;
320 else if (chan == 1)
321 gam = 2.0;
322 else if (chan == 2)
323 gam = 2.1;
324 else
325 gam = 1.0;
326 rv = pow(rv, 1.0/gam);
327 #else
328 rv = icxInvSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], rv,
329 p->in_min[chan], p->in_max[chan]);
330 #endif
331 }
332 return rv;
333 }
334
335 #ifdef NEVER /* Not currently used */
336 /* Inverse Lookup a value though shape curves */
xfit_invshpcurves(xfit * p,double * out,double * in)337 static void xfit_invshpcurves(xfit *p, double *out, double *in) {
338 int e;
339
340 for (e = 0; e < p->di; e++) {
341 double val = in[e];
342 if (p->tcomb & oc_i)
343 val = icxInvSTransFunc(p->v + p->shp_offs[e], p->iluord[e], val,
344 p->in_min[e], p->in_max[e]);
345 out[e] = val;
346 }
347 }
348 #endif /* NEVER */
349
350 /* - - - - - - - - - - - - - - - - - */
351
352 /* Lookup values through the shaper/matrix/shaper model */
xfit_shmatsh(xfit * p,double * out,double * in)353 static void xfit_shmatsh(xfit *p, double *out, double *in) {
354 double tin[MXDI];
355 int e, f;
356
357 for (e = 0; e < p->di; e++)
358 tin[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], in[e],
359 p->in_min[e], p->in_max[e]);
360
361 icxCubeInterp(p->v + p->mat_off, p->fdi, p->di, out, tin);
362
363 if (p->flags & XFIT_OUT_LAB) {
364 for (f = 0; f < p->fdi; f++)
365 out[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], out[f],
366 p->out_min[f], p->out_max[f]);
367 } else {
368 for (f = 0; f < p->fdi; f++)
369 out[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], out[f],
370 p->out_min[f], p->out_max[f]);
371 }
372 }
373
374 /* - - - - - - - - - - - - - - - - - - - - */
375 /* Combined input positioning & shaper transfer curve functions */
376
377 //int db = 0;
378
379 /* Lookup a value though the input positioning and shaper curves */
xfit_inpscurve(xfit * p,double in,int chan)380 static double xfit_inpscurve(xfit *p, double in, int chan) {
381 double rv;
382 /* Just shaper curve */
383 if ((p->tcomb & oc_ip) == oc_i) {
384 rv = icxSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], in,
385 p->in_min[chan], p->in_max[chan]);
386
387 /* shaper and positioning */
388 } else if ((p->tcomb & oc_ip) == oc_ip) {
389 double nin, npind; /* normalized in, normalized position in' */
390 int six; /* Span index */
391 double npind0, npind1; /* normalized position in' span values */
392 double npin0, npin1; /* normalized position in span values */
393 double nsind0, nsind1; /* normalized shaper in' span values */
394 double nsind; /* normalised shaper in' value */
395
396 /* Normalize */
397 nin = (in - p->in_min[chan])/(p->in_max[chan] - p->in_min[chan]);
398
399 //if (db) printf("\n~1 inpscurve: cha %d, input value %f, norm %f\n",chan,in,nin);
400
401 /* Locate the span the input point will be in after positioning lookup */
402 npind = icxTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], nin);
403
404 /* Quantize position space value to grid */
405 six = (int)floor(npind * (p->gres[chan]-1.0));
406 if (six > (p->gres[chan]-2))
407 six = (p->gres[chan]-2);
408
409 /* Compute span position position in' values */
410 npind0 = six / (p->gres[chan]-1.0);
411 npind1 = (six + 1.0) / (p->gres[chan]-1.0);
412
413 //if (db) printf("~1 npind %f, six %d, npind0 %f, npind1 %f\n",npind,six,npind0,npind1);
414
415 /* Compute span in values */
416 npin0 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind0);
417 npin1 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind1);
418
419 //if (db) printf("~1 npin0 %f, npin1 %f\n",npin0,npin1);
420
421 /* Compute shaper space values of in' and spane */
422 #ifdef NEVER
423 nsind = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], nin);
424 nsind0 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin0);
425 nsind1 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin1);
426 #else
427 nsind = xfit_shpcurve(p, nin, chan);
428 nsind0 = xfit_shpcurve(p, npin0, chan);
429 nsind1 = xfit_shpcurve(p, npin1, chan);
430 #endif
431
432 //if (db) printf("~1 nsind %f, nsind0 %f, nsind1 %f\n",nsind,nsind0,nsind1);
433
434 /* Offset and scale shaper in' value to match position span */
435 rv = (nsind - nsind0)/(nsind1 - nsind0) * (npind1 - npind0) + npind0;
436
437 //if (db) printf("~1 scale offset ind %f\n",rv);
438
439 /* de-normalize */
440 rv = rv * (p->in_max[chan] - p->in_min[chan]) + p->in_min[chan];
441 //if (db) printf("~1 returning %d\n",rv);
442
443 /* Just positioning curve */
444 } else if ((p->tcomb & oc_ip) == oc_p) {
445 rv = icxSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], in,
446 p->in_min[chan], p->in_max[chan]);
447 } else {
448 rv = in;
449 }
450 return rv;
451 }
452
453 /* Lookup di values though the input positioning and shaper curves */
xfit_inpscurves(xfit * p,double * out,double * in)454 static void xfit_inpscurves(xfit *p, double *out, double *in) {
455 int e;
456
457 for (e = 0; e < p->di; e++)
458 out[e] = xfit_inpscurve(p, in[e], e);
459 }
460
461 /* Inverse Lookup a value though the input positioning and shaper curves */
xfit_invinpscurve(xfit * p,double in,int chan)462 static double xfit_invinpscurve(xfit *p, double in, int chan) {
463 double rv;
464 /* Just shaper curve */
465 if ((p->tcomb & oc_ip) == oc_i) {
466 rv = icxInvSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], in,
467 p->in_min[chan], p->in_max[chan]);
468
469 /* shaper and positioning */
470 } else if ((p->tcomb & oc_ip) == oc_ip) {
471 double nind, nin; /* normalized in', normalized in */
472 int six; /* Span index */
473 double npind0, npind1; /* normalized position in' span values */
474 double npin0, npin1; /* normalized position in span values */
475 double nsind0, nsind1; /* normalized shaper in' span values */
476 double nsind; /* normalized shaper in' value */
477
478 /* Normalize */
479 nind = (in - p->in_min[chan])/(p->in_max[chan] - p->in_min[chan]);
480
481 //if (db) printf("\n~1 invinpscurve: cha %d, input value %f, norm %f\n",chan,in,nind);
482
483 /* Quantize to grid */
484 six = (int)floor(nind * (p->gres[chan]-1.0));
485 if (six > (p->gres[chan]-2))
486 six = (p->gres[chan]-2);
487
488 /* Compute span in' values */
489 npind0 = six / (p->gres[chan]-1.0);
490 npind1 = (six + 1.0) / (p->gres[chan]-1.0);
491
492 //if (db) printf("~1 six %d, npind0 %f, npind1 %f\n",six,npind0,npind1);
493
494 /* Lookup span in values through position curve */
495 npin0 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind0);
496 npin1 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind1);
497
498 //if (db) printf("~1 npin0 %f, npin1 %f\n",npin0,npin1);
499
500 /* Compute span shaper in' values */
501 #ifdef NEVER
502 nsind0 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin0);
503 nsind1 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin1);
504 #else
505 nsind0 = xfit_shpcurve(p, npin0, chan);
506 nsind1 = xfit_shpcurve(p, npin1, chan);
507 #endif
508
509 /* Offset and scale position in' value to match shaper span */
510 nsind = (nind - npind0)/(npind1 - npind0) * (nsind1 - nsind0) + nsind0;
511
512 //if (db) printf("~1 nsind %f, nsind0 %f, nsind1 %f\n",nsind,nsind0,nsind1);
513
514 /* Invert through shaper curve */
515 #ifdef NEVER
516 nin = icxInvTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], nsind);
517 #else
518 nin = xfit_invshpcurve(p, nsind, chan);
519 #endif
520
521 /* de-normalize */
522 rv = nin * (p->in_max[chan] - p->in_min[chan]) + p->in_min[chan];
523
524 //if (db) printf("\n~1 nin = %f, returning %f\n",nin,rv);
525
526 /* Just positioning curve */
527 } else if ((p->tcomb & oc_ip) == oc_p) {
528 rv = icxInvSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], in,
529 p->in_min[chan], p->in_max[chan]);
530 } else {
531 rv = in;
532 }
533 return rv;
534 }
535
536 #ifdef NEVER
537 /* Check that inverse is working */
_xfit_inpscurve(xfit * p,double in,int chan)538 static double _xfit_inpscurve(xfit *p, double in, int chan) {
539 double inv, rv, iinv;
540
541 inv = in;
542 rv = xfit_inpscurve(p, in, chan);
543 iinv = xfit_invinpscurve(p, rv, chan);
544
545 if (fabs(in - iinv) > 1e-5)
546 warning("xfit_inpscurve check, got %f, should be %f\n",iinv,in);
547
548 return rv;
549 }
550 #define xfit_inpscurve _xfit_inpscurve
551 #endif /* NEVER */
552
553 /* Inverse Lookup di values though the input positioning and shaper curves */
xfit_invinpscurves(xfit * p,double * out,double * in)554 static void xfit_invinpscurves(xfit *p, double *out, double *in) {
555 int e;
556
557 for (e = 0; e < p->di; e++)
558 out[e] = xfit_invinpscurve(p, in[e], e);
559 }
560
561 /* - - - - - - - - - - - - - - - - - - - - */
562 /* Combined output transfer curve functions */
563
564 /* Lookup a value though an output curve */
xfit_outcurve(xfit * p,double in,int chan)565 static double xfit_outcurve(xfit *p, double in, int chan) {
566 double rv;
567 if (p->tcomb & oc_o) {
568 if (p->flags & XFIT_OUT_LAB)
569 rv = icxSTransFunc(p->v + p->out_offs[chan], p->oluord[chan], in,
570 p->out_min[chan], p->out_max[chan]);
571 else
572 rv = icxSTransFuncY2L(p->v + p->out_offs[chan], p->oluord[chan], in,
573 p->out_min[chan], p->out_max[chan]);
574 } else
575 rv = in;
576 return rv;
577 }
578
579 /* Lookup fdi values though the output curves */
xfit_outcurves(xfit * p,double * out,double * in)580 static void xfit_outcurves(xfit *p, double *out, double *in) {
581 int f;
582
583 if (p->flags & XFIT_OUT_LAB) {
584 for (f = 0; f < p->fdi; f++) {
585 double val = in[f];
586 if (p->tcomb & oc_o)
587 val = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], val,
588 p->out_min[f], p->out_max[f]);
589 out[f] = val;
590 }
591 } else {
592 for (f = 0; f < p->fdi; f++) {
593 double val = in[f];
594 if (p->tcomb & oc_o)
595 val = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], val,
596 p->out_min[f], p->out_max[f]);
597 out[f] = val;
598 }
599 }
600 }
601
602 /* Inverse Lookup a value though an output curve */
xfit_invoutcurve(xfit * p,double in,int chan)603 static double xfit_invoutcurve(xfit *p, double in, int chan) {
604 double rv;
605 if (p->tcomb & oc_o) {
606 if (p->flags & XFIT_OUT_LAB)
607 rv = icxInvSTransFunc(p->v + p->out_offs[chan], p->oluord[chan], in,
608 p->out_min[chan], p->out_max[chan]);
609 else
610 rv = icxInvSTransFuncY2L(p->v + p->out_offs[chan], p->oluord[chan], in,
611 p->out_min[chan], p->out_max[chan]);
612 } else
613 rv = in;
614 return rv;
615 }
616
617 /* Inverse Lookup fdi values though output curves */
xfit_invoutcurves(xfit * p,double * out,double * in)618 static void xfit_invoutcurves(xfit *p, double *out, double *in) {
619 int f;
620
621 if (p->flags & XFIT_OUT_LAB) {
622 for (f = 0; f < p->fdi; f++) {
623 double val = in[f];
624 if (p->tcomb & oc_o)
625 val = icxInvSTransFunc(p->v + p->out_offs[f], p->oluord[f], val,
626 p->out_min[f], p->out_max[f]);
627 out[f] = val;
628 }
629 } else {
630 for (f = 0; f < p->fdi; f++) {
631 double val = in[f];
632 if (p->tcomb & oc_o)
633 val = icxInvSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], val,
634 p->out_min[f], p->out_max[f]);
635 out[f] = val;
636 }
637 }
638 }
639
640 /* - - - - - - - - - */
641
642 /* Convert an output value from absolute */
643 /* to relative using the current white point. */
xfit_abs_to_rel(xfit * p,double * out,double * in)644 static void xfit_abs_to_rel(xfit *p, double *out, double *in) {
645 if (p->flags & XFIT_OUT_WP_REL) {
646 if (p->flags & XFIT_OUT_LAB) {
647 icmLab2XYZ(&icmD50, out, in);
648 icmMulBy3x3(out, p->fromAbs, out);
649 icmXYZ2Lab(&icmD50, out, out);
650 } else {
651 icmMulBy3x3(out, p->fromAbs, in);
652 }
653 } else {
654 out[0] = in[0];
655 out[1] = in[1];
656 out[2] = in[2];
657 }
658 }
659
660 /* Convert an XYZ output value from absolute */
661 /* to cLut relative using the current white point. */
xfit_XYZ_abs_to_rel(xfit * p,double * out,double * in)662 static void xfit_XYZ_abs_to_rel(xfit *p, double *out, double *in) {
663 if (p->flags & XFIT_OUT_WP_REL) {
664 if (p->flags & XFIT_OUT_LAB) {
665 icmMulBy3x3(out, p->fromAbs, in);
666 icmXYZ2Lab(&icmD50, out, out);
667 } else {
668 icmMulBy3x3(out, p->fromAbs, in);
669 }
670 } else {
671 if (p->flags & XFIT_OUT_LAB) {
672 icmXYZ2Lab(&icmD50, out, in);
673 } else {
674 out[0] = in[0];
675 out[1] = in[1];
676 out[2] = in[2];
677 }
678 }
679 }
680
681 /* - - - - - - - - - */
682
683 /* return a weighting for the magnitude of the in and out */
684 /* shaping parameters squared. This is to reduce unconstrained "wiggles" */
shapmag(xfit * p)685 static double shapmag(
686 xfit *p /* Base of optimisation structure */
687 ) {
688 double tt, w;
689 double *b; /* Base of parameters for this section */
690 int di = p->di;
691 int fdi = p->fdi;
692 int e, f, k;
693 double iparam = 0.0;
694 double oparam = 0.0;
695 double dd;
696
697 if (p->opt_msk & oc_i) {
698 dd = SHAPE_WEIGHT/(double)(di);
699 b = p->v + p->shp_off;
700 for (e = 0; e < di; e++) {
701 for (k = 0; k < p->iluord[e]; k++) {
702 if (k <= 1) {
703 w = SHAPE_HW01;
704 } else if (k <= SHAPE_HBREAK) {
705 double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
706 w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
707 w *= p->shp_smooth[e];
708 } else {
709 w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
710 w *= p->shp_smooth[e];
711 }
712 tt = *b++;
713 tt *= tt; /* Squared */
714 iparam += w * tt;
715 }
716 }
717 iparam *= dd;
718 }
719
720 if (p->opt_msk & oc_o) {
721 dd = SHAPE_WEIGHT/(double)(fdi);
722 b = p->v + p->out_off;
723 for (f = 0; f < fdi; f++) {
724 for (k = 0; k < p->oluord[f]; k++) {
725 if (k <= 1) {
726 w = SHAPE_HW01;
727 } else if (k <= SHAPE_HBREAK) {
728 double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
729 w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
730 w *= p->out_smooth[f];
731 } else {
732 w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
733 w *= p->out_smooth[f];
734 }
735 tt = *b++;
736 tt *= tt; /* Squared */
737 oparam += w * tt;
738 }
739 }
740 oparam *= dd;
741 }
742 return iparam + oparam;
743 }
744
745 /* return a weighting for the magnitude of the in and out */
746 /* shaping parameters. This is to reduce unconstrained "wiggles" */
747 /* Also sum the partial derivative for the parameters involved */
dshapmag(xfit * p,double * dav)748 static double dshapmag(
749 xfit *p, /* Base of optimisation structure */
750 double *dav /* Sum del's */
751 ) {
752 double tt, w;
753 double *b, *c; /* Base of parameters for this section */
754 int di = p->di;
755 int fdi = p->fdi;
756 int e, f, k;
757 double iparam = 0.0;
758 double oparam = 0.0;
759 double dd;
760
761 if (p->opt_msk & oc_i) {
762 dd = SHAPE_WEIGHT/(double)(di);
763 b = p->v + p->shp_off;
764 c = dav + p->shp_off;
765 for (e = 0; e < di; e++) {
766 for (k = 0; k < p->iluord[e]; k++) {
767 if (k <= 1) {
768 w = SHAPE_HW01;
769 } else if (k <= SHAPE_HBREAK) {
770 double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
771 w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
772 w *= p->shp_smooth[e];
773 } else {
774 w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
775 w *= p->shp_smooth[e];
776 }
777 tt = *b++;
778 *c++ += 2.0 * dd * w * tt;
779 tt *= tt; /* Squared */
780 iparam += w * tt;
781
782 }
783 }
784 iparam *= dd;
785 }
786
787 if (p->opt_msk & oc_o) {
788 dd = SHAPE_WEIGHT/(double)(fdi);
789 b = p->v + p->out_off;
790 c = dav + p->out_off;
791 for (f = 0; f < fdi; f++) {
792 for (k = 0; k < p->oluord[f]; k++) {
793 if (k <= 1) {
794 w = SHAPE_HW01;
795 } else if (k <= SHAPE_HBREAK) {
796 double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
797 w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
798 w *= p->out_smooth[f];
799 } else {
800 w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
801 w *= p->out_smooth[f];
802 }
803 tt = *b++;
804 *c++ += 2.0 * dd * w * tt;
805 tt *= tt; /* Squared */
806 oparam += w * tt;
807 }
808 }
809 oparam *= dd;
810 }
811 return iparam + oparam;
812 }
813
814 /* Scale the shaper derivatives */
dshapscale(xfit * p,double * dav,double scale)815 static void dshapscale(
816 xfit *p, /* Base of optimisation structure */
817 double *dav, /* del's */
818 double scale /* Scale factor */
819 ) {
820 double tt, w;
821 double *b, *c; /* Base of parameters for this section */
822 int di = p->di;
823 int fdi = p->fdi;
824 int e, f, k;
825
826 if (p->opt_msk & oc_i) {
827 c = dav + p->shp_off;
828 for (e = 0; e < di; e++) {
829 for (k = 0; k < p->iluord[e]; k++) {
830 *c++ *= scale;
831 }
832 }
833 }
834
835 if (p->opt_msk & oc_o) {
836 c = dav + p->out_off;
837 for (f = 0; f < fdi; f++) {
838 for (k = 0; k < p->oluord[f]; k++) {
839 *c++ *= scale;
840 }
841 }
842 }
843 }
844
845 /* Progress function */
xfitprog(void * pdata,int perc)846 static void xfitprog(void *pdata, int perc) {
847 xfit *p = (xfit *)pdata;
848
849 if (p->verb) {
850 printf("%c% 3d%%",cr_char,perc);
851 if (perc == 100)
852 printf("\n");
853 fflush(stdout);
854 }
855 }
856
857
858 int xfitfunc_trace = 1;
859
860 /* Shaper+Matrix optimisation function handed to powell() */
861 /* We simply minimize the total delta E squared, consistent with smoothness */
xfitfunc(void * edata,double * v)862 static double xfitfunc(void *edata, double *v) {
863 xfit *p = (xfit *)edata;
864 double tw = 0.0; /* Total weight */
865 double ev = 0.0, rv, smv;
866 double tin[MXDI], out[MXDO];
867 int di = p->di;
868 int fdi = p->fdi;
869 int i, e, f;
870
871 /* Copy the parameters being optimised into xfit structure */
872
873 /* Special case - a single shaper curve. The first sm_iluord params */
874 /* are the common curve parameters, and the remainder are the matrix onwards */
875 if (p->opt_ssch) {
876
877 for (e = 0; e < di; e++) { /* Duplicate and extend to per channel curve params */
878 for (i = 0; i < p->sm_iluord; i++)
879 p->v[p->shp_offs[e] + i] = v[i];
880 for (; i < p->iluord[e]; i++)
881 p->v[p->shp_offs[e] + i] = 0.0;
882 }
883 for (i = p->sm_iluord; i < p->opt_cnt; i++)
884 p->v[p->mat_off + i - p->sm_iluord] = v[i];
885 } else {
886 for (i = 0; i < p->opt_cnt; i++) {
887 //printf("~1 param %d = %f\n",i,v[i]);
888 p->v[p->opt_off + i] = v[i];
889 }
890 }
891
892 /* For all our data points */
893 for (i = 0; i < p->nodp; i++) {
894 double del;
895
896 /* Apply input shaper channel curves */
897 for (e = 0; e < di; e++)
898 tin[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], p->rpoints[i].p[e],
899 p->in_min[e], p->in_max[e]);
900
901 /* Apply matrix cube interpolation */
902 icxCubeInterp(p->v + p->mat_off, fdi, di, out, tin);
903
904 /* Apply output channel curves */
905 for (f = 0; f < fdi; f++) {
906 if (p->flags & XFIT_OUT_LAB) {
907 out[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], out[f],
908 p->out_min[f], p->out_max[f]);
909 } else {
910 out[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], out[f],
911 p->out_min[f], p->out_max[f]);
912 }
913 }
914
915 /* Evaluate the error squared */
916 if (p->flags & XFIT_FM_INPUT) {
917 double pp[MXDI];
918 for (e = 0; e < di; e++)
919 pp[e] = p->rpoints[i].p[e];
920 for (f = 0; f < fdi; f++) {
921 double t1 = p->rpoints[i].v[f] - out[f]; /* Error in output */
922
923 /* Create input point offset by equivalent delta to output point */
924 /* error, in proportion to the partial derivatives for that output. */
925 for (e = 0; e < di; e++)
926 pp[e] += t1 * p->piv[i].ide[f][e];
927 }
928 del = p->to_de2(p->cntx2, pp, p->rpoints[i].p);
929 } else {
930 del = p->to_de2(p->cntx2, out, p->rpoints[i].v);
931 }
932 if (CURVEPOW > 1.0)
933 del = pow(del, CURVEPOW);
934 tw += p->rpoints[i].w;
935 ev += p->rpoints[i].w * del;
936 }
937
938 /* Normalise error to be an average delta E squared */
939 ev /= tw;
940
941 /* Sum with shaper parameters squared, to */
942 /* minimise unsconstrained "wiggles" */
943 smv = shapmag(p);
944 if (CURVEPOW > 1.0)
945 smv = pow(smv, CURVEPOW);
946 rv = ev + smv;
947
948 #ifdef DEBUG_PROGRESS
949 if (xfitfunc_trace)
950 fprintf(stdout,"~1(sm %f, ev %f)xfitfunc returning %f\n",smv,ev,rv);
951 #endif
952
953 return rv;
954 }
955
956 /* Shaper+Matrix optimisation function with partial derivatives, */
957 /* handed to conjgrad() */
dxfitfunc(void * edata,double * dv,double * v)958 static double dxfitfunc(void *edata, double *dv, double *v) {
959 xfit *p = (xfit *)edata;
960 double tw = 0.0; /* Total weight */
961 double ev = 0.0, rv, smv;
962 double tin[MXDI], out[MXDO];
963
964 double dav[MXPARMS]; /* Overall del due to del param vals */
965 double sdav[MXPARMS]; /* Overall del due to del smooth param vals */
966
967 double dtin_iv[MXDI * MXLUORD]; /* Del in itrans out due to del itrans param vals */
968 double dmato_mv[1 << MXDI]; /* Del in mat out due to del in matrix param vals */
969 double dmato_tin[MXDO * MXDI]; /* Del in mat out due to del in matrix input values */
970 double dout_ov[MXDO * MXLUORD]; /* Del in otrans out due to del in otrans param values */
971 double dout_mato[MXDO]; /* Del in otrans out due to del in otrans input values */
972
973 double dout_de[2][MXDIDO]; /* Del in DE due to two output values */
974
975 int di = p->di;
976 int fdi = p->fdi;
977 int i, jj, k, e, ee, f, ff;
978
979 /* Copy the parameters being optimised into xfit structure */
980
981 /* Special case - a single shaper curve. The first sm_iluord params */
982 /* are the common curve parameters, and the remainder are the matrix onwards */
983 if (p->opt_ssch) {
984 for (e = 0; e < di; e++) { /* Duplicate and extend to per channel curve params */
985 for (i = 0; i < p->sm_iluord; i++)
986 p->v[p->shp_offs[e] + i] = v[i];
987 for (; i < p->iluord[e]; i++)
988 p->v[p->shp_offs[e] + i] = 0.0;
989 }
990 for (i = p->sm_iluord; i < p->opt_cnt; i++)
991 p->v[p->mat_off + i - p->sm_iluord] = v[i];
992
993 } else {
994 for (i = 0; i < p->opt_cnt; i++) {
995 p->v[p->opt_off + i] = v[i];
996 }
997 }
998
999 /* Zero the accumulated partial derivatives */
1000 /* We compute deriv for all parameters (not just current optimised) */
1001 for (i = 0; i < p->tot_cnt; i++)
1002 dav[i] = 0.0;
1003
1004 /* For all our data points */
1005 for (i = 0; i < p->nodp; i++) {
1006 double del;
1007
1008 /* Apply input channel curves */
1009 for (e = 0; e < di; e++)
1010 tin[e] = icxdpSTransFunc(p->v + p->shp_offs[e], &dtin_iv[p->shp_offs[e] - p->shp_off],
1011 p->iluord[e], p->rpoints[i].p[e], p->in_min[e], p->in_max[e]);
1012
1013 /* Apply matrix cube interpolation */
1014 icxdpdiCubeInterp(p->v + p->mat_off, dmato_mv, dmato_tin, fdi, di, out, tin);
1015
1016 /* Apply output channel curves */
1017 for (f = 0; f < fdi; f++) {
1018 if (p->flags & XFIT_OUT_LAB)
1019 out[f] = icxdpdiSTransFunc(p->v + p->out_offs[f],
1020 &dout_ov[p->out_offs[f] - p->out_off], &dout_mato[f],
1021 p->oluord[f], out[f], p->out_min[f], p->out_max[f]);
1022 else
1023 out[f] = icxdpdiSTransFuncY2L(p->v + p->out_offs[f],
1024 &dout_ov[p->out_offs[f] - p->out_off], &dout_mato[f],
1025 p->oluord[f], out[f], p->out_min[f], p->out_max[f]);
1026 }
1027
1028 /* Convert to Delta E and compute pde's into dout_de squared */
1029 if (p->flags & XFIT_FM_INPUT) {
1030 double tdout_de[2][MXDIDO];
1031 double pp[MXDI];
1032 for (e = 0; e < di; e++)
1033 pp[e] = p->rpoints[i].p[e];
1034 for (f = 0; f < fdi; f++) {
1035 double t1 = p->rpoints[i].v[f] - out[f]; /* Error in output */
1036
1037 /* Create input point offset by equivalent delta to output point */
1038 /* error, in proportion to the partial derivatives for that output. */
1039 for (e = 0; e < di; e++)
1040 pp[e] += t1 * p->piv[i].ide[f][e];
1041 }
1042 del = p->to_dde2(p->cntx2, tdout_de, pp, p->rpoints[i].p);
1043 if (CURVEPOW > 1.0) {
1044 double dadj;
1045 dadj = CURVEPOW * pow(del, CURVEPOW - 1.0); /* Adjust derivative accordingly */
1046 del = pow(del, CURVEPOW);
1047 for (e = 0; e < di; e++)
1048 tdout_de[0][e] *= dadj;
1049 }
1050
1051 /* Compute partial derivative */
1052 for (e = 0; e < di; e++) {
1053 dout_de[0][e] = 0.0;
1054 for (f = 0; f < fdi; f++) {
1055 dout_de[0][e] += tdout_de[0][e] * p->piv[i].ide[f][e];
1056 }
1057 }
1058 } else {
1059 del = p->to_dde2(p->cntx2, dout_de, out, p->rpoints[i].v);
1060 if (CURVEPOW > 1.0) {
1061 double dadj;
1062 dadj = CURVEPOW * pow(del, CURVEPOW - 1.0); /* Adjust derivative accordingly */
1063 del = pow(del, CURVEPOW);
1064 for (f = 0; f < fdi; f++)
1065 dout_de[0][f] *= dadj;
1066 }
1067 }
1068
1069 /* Accumulate total weighted delta E squared */
1070 tw += p->rpoints[i].w;
1071 ev += p->rpoints[i].w * del;
1072
1073 /* Compute and accumulate partial difference values for each parameter value */
1074 if (p->opt_msk & oc_i) {
1075 /* Input transfer parameters */
1076 for (ee = 0; ee < di; ee++) { /* Parameter input chanel */
1077 for (k = 0; k < p->iluord[ee]; k++) { /* Param within channel */
1078 double vv = 0.0;
1079 jj = p->shp_offs[ee] - p->shp_off + k; /* Overall input trans param */
1080
1081 // for (ff = 0; ff < 3; ff++) { /* Lab channels */
1082 for (ff = 0; ff < fdi; ff++) { /* Output channels */
1083 vv += dout_de[0][ff] * dout_mato[ff]
1084 * dmato_tin[ff * di + ee] * dtin_iv[jj];
1085 }
1086 dav[p->shp_off + jj] += p->rpoints[i].w * vv;
1087 }
1088 }
1089 }
1090
1091 if (p->opt_msk & oc_m) {
1092 /* Matrix parameters */
1093 for (ff = 0; ff < fdi; ff++) { /* Parameter output chanel */
1094 for (ee = 0; ee < (1 << di); ee++) { /* Matrix input combination chanel */
1095 double vv = 0.0;
1096 jj = ff * (1 << di) + ee; /* Matrix Parameter index */
1097
1098 vv += dout_de[0][ff] * dout_mato[ff] * dmato_mv[ee];
1099 dav[p->mat_off + jj] += p->rpoints[i].w * vv;
1100 }
1101 }
1102 }
1103
1104 if (p->opt_msk & oc_o) {
1105 /* Output transfer parameters */
1106 for (ff = 0; ff < fdi; ff++) { /* Parameter output chanel */
1107 for (k = 0; k < p->oluord[ff]; k++) { /* Param within channel */
1108 double vv = 0.0;
1109 jj = p->out_offs[ff] - p->out_off + k; /* Overall output trans param */
1110
1111 vv += dout_de[0][ff] * dout_ov[jj];
1112 dav[p->out_off + jj] += p->rpoints[i].w * vv;
1113 }
1114 }
1115 }
1116 }
1117
1118 /* Normalise error to be an average delta E squared */
1119 ev /= tw;
1120 for (i = 0; i < p->tot_cnt; i++) {
1121 dav[i] /= tw;
1122 sdav[i] = 0.0;
1123 }
1124
1125 /* Sum with shaper parameters squared, to */
1126 /* minimise unsconstrained "wiggles" */
1127 /* Compute partial derivative wrt those parameters too */
1128 smv = dshapmag(p, sdav);
1129 if (CURVEPOW > 1.0) {
1130 double dadj;
1131 dadj = CURVEPOW * pow(smv, CURVEPOW - 1.0); /* Adjust derivative accordingly */
1132 smv = pow(smv, CURVEPOW);
1133 dshapscale(p, sdav, dadj); /* Scale the partial derivatives */
1134 }
1135 rv = ev + smv;
1136
1137 /* Sum the del for parameters being optimised and copy to return array */
1138 if (p->opt_ssch) {
1139 for (i = 0; i < p->sm_iluord; i++)
1140 dv[i] = 0.0;
1141 for (e = 0; e < di; e++) { /* Combine per channel curve de's */
1142 for (i = 0; i < p->sm_iluord; i++)
1143 dv[i] += dav[p->shp_offs[e] + i] + sdav[p->shp_offs[e] + i];
1144 }
1145 for (i = p->sm_iluord; i < p->opt_cnt; i++) /* matrix and rest de's */
1146 dv[i] = dav[p->mat_off + i - p->sm_iluord] + sdav[p->mat_off + i - p->sm_iluord];
1147
1148 } else {
1149 for (i = 0; i < p->opt_cnt; i++)
1150 dv[i] = dav[p->opt_off + i] + sdav[p->opt_off + i];
1151 }
1152
1153 #ifdef DEBUG_PROGRESS
1154 fprintf(stdout,"~1(sm %f, ev %f)dxfitfunc returning %f\n",smv,ev,rv);
1155 #endif
1156
1157 return rv;
1158 }
1159
1160 #ifdef NEVER
1161 /* Check partial derivative function within xfitfunc() [Intensive check] */
1162
_xfitfunc(void * edata,double * v)1163 static double _xfitfunc(void *edata, double *v) {
1164 xfit *p = (xfit *)edata;
1165 int i;
1166 double dv[MXPARMS];
1167 double rv, drv;
1168 double trv;
1169 int verb;
1170
1171 rv = xfitfunc(edata, v);
1172 verb = p->verb;
1173 p->verb = 0;
1174 drv = dxfitfunc(edata, dv, v);
1175 p->verb = verb;
1176
1177 if (fabs(rv - drv) > 1e-6)
1178 printf("######## RV MISMATCH is %f should be %f ########\n",rv,drv);
1179
1180 /* Check each parameter delta */
1181 xfitfunc_trace = 0;
1182 for (i = 0; i < p->opt_cnt; i++) {
1183 double del;
1184
1185 v[i] += 1e-7;
1186 trv = xfitfunc(edata, v);
1187 v[i] -= 1e-7;
1188
1189 /* Check that del is correct */
1190 del = (trv - rv)/1e-7;
1191 if (fabs(dv[i] - del) > 0.04) {
1192 //printf("~1 del = %f from (trv %f - rv %f)/0.1\n",del,trv,rv);
1193 printf("######## EXCESSIVE at v[%d] is %f should be %f ########\n",i,dv[i],del);
1194 }
1195 }
1196 xfitfunc_trace = 1;
1197 return rv;
1198 }
1199
1200 #define xfitfunc _xfitfunc
1201 #endif /* NEVER */
1202
1203 #ifdef NEVER
1204 /* Check partial derivative function within dxfitfunc() [Less intensive check] */
1205
_dxfitfunc(void * edata,double * dv,double * v)1206 static double _dxfitfunc(void *edata, double *dv, double *v) {
1207 xfit *p = (xfit *)edata;
1208 int i;
1209 double rv, drv;
1210 double trv;
1211 int verb;
1212 int exec = 0;
1213
1214 rv = xfitfunc(edata, v);
1215 verb = p->verb;
1216 p->verb = 0;
1217 drv = dxfitfunc(edata, dv, v);
1218 p->verb = verb;
1219
1220 if (fabs(rv - drv) > 1e-6)
1221 printf("######## RV MISMATCH is %f should be %f ########\n",rv,drv);
1222
1223 /* Check each parameter delta */
1224 xfitfunc_trace = 0;
1225 for (i = 0; i < p->opt_cnt; i++) {
1226 double del;
1227
1228 v[i] += 1e-7;
1229 trv = xfitfunc(edata, v);
1230 v[i] -= 1e-7;
1231
1232 /* Check that del is correct */
1233 del = (trv - rv)/1e-7;
1234 if (fabs(dv[i] - del) > 0.04) {
1235 //printf("~1 del = %f from (trv %f - rv %f)/0.1\n",del,trv,rv);
1236 printf("######## EXCESSIVE at v[%d] is %f should be %f ########\n",i,dv[i],del);
1237 exec = 1;
1238 }
1239 }
1240 #ifdef NEVER
1241 if (exec) {
1242 printf("~1 parameters are:\n");
1243 for (i = 0; i < p->opt_cnt; i++)
1244 printf("p->wv[%d] = %f;\n",i,v[i]);
1245 exit(1);
1246 }
1247 #endif
1248 xfitfunc_trace = 1;
1249 return rv;
1250 }
1251
1252 #define dxfitfunc _dxfitfunc
1253 #endif /* NEVER */
1254
1255 /* - - - - - - - - - */
1256
1257 /* Output curve symetry optimisation function handed to powell() */
1258 /* Just the order 0 value will be adjusted */
symoptfunc(void * edata,double * v)1259 static double symoptfunc(void *edata, double *v) {
1260 xfit *p = (xfit *)edata;
1261 double out[1], in[1] = { 0.0 };
1262 int ch = p->opt_ch; /* Output channel being adjusted for symetry */
1263 double rv;
1264
1265 /* Copy the parameter being tested back into xfit */
1266 p->v[p->out_offs[ch]] = v[0];
1267 if (p->flags & XFIT_OUT_LAB)
1268 *out = icxSTransFunc(p->v + p->out_offs[ch], p->oluord[ch], *in,
1269 p->out_min[ch], p->out_max[ch]);
1270 else
1271 *out = icxSTransFuncY2L(p->v + p->out_offs[ch], p->oluord[ch], *in,
1272 p->out_min[ch], p->out_max[ch]);
1273
1274 rv = out[0] * out[0];
1275
1276 #ifdef DEBUG_PROGRESS
1277 printf("~1symoptfunc returning %f\n",rv);
1278 #endif
1279 return rv;
1280 }
1281
1282 /* - - - - - - - - - */
1283
1284 /* Set up for an optimisation run: */
1285 /* Figure out parameters being optimised, */
1286 /* copy them to start values, */
1287 /* init and scale the search radius */
setup_xfit(xfit * p,double * wv,double * sa,double transrad,double matrad)1288 static void setup_xfit(
1289 xfit *p,
1290 double *wv, /* Return parameters to hand to optimiser */
1291 double *sa, /* Return search radius to hand to optimiser */
1292 double transrad,/* Nominal transfer curve radius, 0.0 - 3.0 */
1293 double matrad /* Nominal matrix radius, 0.0 - 1.0 */
1294 ) {
1295 int i;
1296 p->opt_off = -1;
1297 p->opt_cnt = 0;
1298
1299 if (p->opt_msk & oc_i) {
1300
1301 if (p->opt_ssch) { /* Special case - should only be used first, */
1302 /* Fitting a sigle common input shaper curve. */
1303 if (p->opt_off < 0)
1304 p->opt_off = p->mat_off - p->sm_iluord; /* Shouldn't be used... */
1305 p->opt_cnt += p->sm_iluord;
1306
1307 for (i = 0; i < p->sm_iluord; i++) {
1308 *wv++ = 0.0;
1309 *sa++ = transrad;
1310 }
1311
1312 } else { /* Initial or continuing fitting of all the curves */
1313 if (p->opt_off < 0)
1314 p->opt_off = p->shp_off;
1315 p->opt_cnt += p->shp_cnt;
1316
1317 for (i = 0; i < p->shp_cnt; i++) {
1318 *wv++ = p->v[p->shp_off + i];
1319 *sa++ = transrad;
1320 }
1321 }
1322 }
1323 if (p->opt_msk & oc_m) {
1324 if (p->opt_off < 0)
1325 p->opt_off = p->mat_off;
1326 p->opt_cnt += p->mat_cnt;
1327
1328 for (i = 0; i < p->mat_cnt; i++) {
1329 *wv++ = p->v[p->mat_off + i];
1330 *sa++ = matrad;
1331 }
1332 }
1333 if (p->opt_msk & oc_o) {
1334 if (p->opt_off < 0)
1335 p->opt_off = p->out_off;
1336 p->opt_cnt += p->out_cnt;
1337
1338 for (i = 0; i < p->out_cnt; i++) {
1339 *wv++ = p->v[p->out_off + i];
1340 *sa++ = transrad;
1341 }
1342 }
1343 if (p->opt_cnt > MXPARMS)
1344 error("setup_xfit: asert, %d exceeded MXPARMS %d",p->opt_cnt,MXPARMS);
1345
1346 #ifdef DEBUG
1347 printf("setup_xfit() got opt_msk 0x%x, opt_off %d, opt_cnt = %d\n",
1348 p->opt_msk,p->opt_off,p->opt_cnt);
1349 #endif /* DEBUG */
1350 }
1351
1352 #ifdef DEBUG
1353 /* Diagnostic */
dump_xfit(xfit * p)1354 static void dump_xfit(
1355 xfit *p
1356 ) {
1357 int i, e, f;
1358 double *b; /* Base of parameters for this section */
1359 int di, fdi;
1360 di = p->di;
1361 fdi = p->fdi;
1362
1363 /* Input positioning curve */
1364 b = p->v + p->pos_off;
1365 for (e = 0; e < di; b += p->iluord[e], e++) {
1366 printf("pos %d = ",e);
1367 for (i = 0; i < p->iluord[e]; i++)
1368 printf("%f ",b[i]);
1369 printf("\n");
1370 }
1371
1372 /* Input shaper curve */
1373 b = p->v + p->shp_off;
1374 for (e = 0; e < di; b += p->iluord[e], e++) {
1375 printf("shp %d = ",e);
1376 for (i = 0; i < p->iluord[e]; i++)
1377 printf("%f ",b[i]);
1378 printf("\n");
1379 }
1380
1381 /* Matrix */
1382 b = p->v + p->mat_off;
1383 for (e = 0; e < (1 << di); e++) {
1384 printf("mx %d = ",e);
1385 for (f = 0; f < fdi; f++)
1386 printf("%f ",*b++);
1387 printf("\n");
1388 }
1389
1390 /* Output curve */
1391 b = p->v + p->out_off;
1392 for (f = 0; f < fdi; b += p->oluord[f], f++) {
1393 printf("out %d = ",f);
1394 for (i = 0; i < p->oluord[f]; i++)
1395 printf("%f ",b[i]);
1396 printf("\n");
1397 }
1398 }
1399 #endif /* DEBUG */
1400
1401 /* - - - - - - - - - */
1402
1403 /* Setup the pseudo inverse information for each test point, */
1404 /* using the current model. */
setup_piv(xfit * p)1405 static void setup_piv(xfit *p) {
1406 int di = p->di;
1407 int fdi = p->fdi;
1408 int i, e, f;
1409
1410 /* Estimate in -> out partial derivatives */
1411 for (i = 0; i < p->nodp; i++) {
1412 double pd[MXDO][MXDI];
1413 double pp[MXDI];
1414 double vv[MXDIDO];
1415
1416 /* Estimate in -> out partial derivatives */
1417 for (e = 0; e < di; e++)
1418 pp[e] = p->ipoints[i].p[e];
1419
1420 /* Apply input shaper channel curves */
1421 for (e = 0; e < di; e++)
1422 vv[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], pp[e],
1423 p->in_min[e], p->in_max[e]);
1424
1425 /* Apply matrix cube interpolation */
1426 icxCubeInterp(p->v + p->mat_off, fdi, di, vv, vv);
1427
1428 /* Apply output channel curves */
1429 for (f = 0; f < fdi; f++) {
1430 if (p->flags & XFIT_OUT_LAB)
1431 vv[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], vv[f],
1432 p->out_min[f], p->out_max[f]);
1433 else
1434 vv[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], vv[f],
1435 p->out_min[f], p->out_max[f]);
1436 }
1437
1438 for (e = 0; e < di; e++) {
1439 double tt[MXDIDO];
1440
1441 pp[e] += 1e-4;
1442
1443 /* Apply input shaper channel curves */
1444 for (e = 0; e < di; e++)
1445 tt[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], pp[e],
1446 p->in_min[e], p->in_max[e]);
1447
1448 /* Apply matrix cube interpolation */
1449 icxCubeInterp(p->v + p->mat_off, fdi, di, tt, tt);
1450
1451 /* Apply output channel curves */
1452 for (f = 0; f < fdi; f++) {
1453 if (p->flags & XFIT_OUT_LAB)
1454 tt[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], tt[f],
1455 p->out_min[f], p->out_max[f]);
1456 else
1457 tt[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], tt[f],
1458 p->out_min[f], p->out_max[f]);
1459 }
1460
1461 for (f = 0; f < p->fdi; f++)
1462 pd[f][e] = (tt[f] - vv[f])/1e-4;
1463
1464 pp[e] -= 1e-4;
1465 }
1466
1467 /* Compute a psudo inverse matrix to map rout delta E to */
1468 /* in delta E in proportion to the pd magnitude. */
1469 for (f = 0; f < fdi; f++) {
1470 double ss = 0.0;
1471
1472 for (e = 0; e < di; e++) /* Sum of pd's ^4 */
1473 ss += pd[f][e] * pd[f][e] * pd[f][e] * pd[f][e];
1474 ss = sqrt(ss);
1475 if (ss > 1e-8) {
1476 for (e = 0; e < di; e++)
1477 p->piv[i].ide[f][e] = pd[f][e]/ss;
1478 } else { /* Hmm. */
1479 for (e = 0; e < di; e++)
1480 p->piv[i].ide[f][e] = 0.0;
1481 }
1482 }
1483 }
1484 }
1485
1486 /* - - - - - - - - - */
1487
1488 /* Function to pass to rspl to re-set output values, */
1489 /* to account for skeleton model offset. */
1490 static void
skm_rspl_out(void * pp,double * out,double * in)1491 skm_rspl_out(
1492 void *pp, /* relativectx structure */
1493 double *out, /* output value */
1494 double *in /* input value */
1495 ) {
1496 xfit *p = (xfit *)pp;
1497 int f, fdi = p->fdi;
1498 double inval[MXDI];
1499 double skval[MXDO];
1500
1501 /* Look up the skeleton value for this grid point */
1502 xfit_invinpscurves(p, inval, in); /* Back to input values */
1503 p->skm->lookup(p->skm, skval, inval); /* Skm */
1504 xfit_abs_to_rel(p, skval, skval);
1505 xfit_invoutcurves(p, skval, skval);
1506
1507 for (f = 0; f < fdi; f++)
1508 out[f] += skval[f]; /* Add it back */
1509 }
1510
1511 /* Weak function rspl callback (not used) */
skm_weak(void * cbntx,double * out,double * in)1512 void skm_weak(void *cbntx, double *out, double *in) {
1513 xfit *p = (xfit *)cbntx;
1514
1515 #ifndef NEVER
1516 int f, fdi = p->fdi;
1517
1518 for (f = 0; f < fdi; f++)
1519 out[f] = 0.0; /* Deviation from skeleton should tend to zero */
1520
1521 #else /* Skeleton as weak atractor */
1522 int f, fdi = p->fdi;
1523 double inval[MXDI];
1524
1525 /* Look up the skeleton value for this grid point */
1526 xfit_invinpscurves(p, inval, in); /* Back to input values */
1527 p->skm->lookup(p->skm, out, inval); /* Skm */
1528 xfit_abs_to_rel(p, out, out);
1529 xfit_invoutcurves(p, out, out);
1530
1531 #endif
1532 }
1533
1534 /* - - - - - - - - - */
1535
1536 /* Function to pass to rspl to re-set output values, */
1537 /* to make them relative to the white and black points */
1538 static void
conv_rspl_out(void * pp,double * out,double * in)1539 conv_rspl_out(
1540 void *pp, /* relativectx structure */
1541 double *out, /* output value */
1542 double *in /* input value */
1543 ) {
1544 xfit *p = (xfit *)pp;
1545 double tt[3];
1546
1547 /* Convert the clut values to output values */
1548 xfit_outcurves(p, tt, out);
1549
1550 if (p->flags & XFIT_OUT_LAB) {
1551 icmLab2XYZ(&icmD50, tt, tt);
1552 icmMulBy3x3(out, p->cmat, tt);
1553 icmXYZ2Lab(&icmD50, out, out);
1554
1555 } else { /* We are all in XYZ */
1556 icmMulBy3x3(out, p->cmat, tt);
1557 }
1558
1559 /* And then convert them back to clut values */
1560 xfit_invoutcurves(p, out, out);
1561 }
1562
1563 /* Function to pass to rspl to re-set output values, */
1564 /* to clip any with Y over 1.0 to D50 */
1565 static void
clip_rspl_out(void * pp,double * out,double * in)1566 clip_rspl_out(
1567 void *pp, /* relativectx structure */
1568 double *out, /* output value */
1569 double *in /* input value */
1570 ) {
1571 xfit *p = (xfit *)pp;
1572 double tt[3];
1573
1574 /* Convert the clut values to output values */
1575 xfit_outcurves(p, tt, out);
1576
1577 if (p->flags & XFIT_OUT_LAB) {
1578 if (tt[0] > 100.0)
1579 icmCpy3(out, p->cmat[0]);
1580 } else {
1581 if (tt[1] > 1.0)
1582 icmCpy3(out, p->cmat[0]);
1583 }
1584 }
1585
1586 //#ifdef SPECIAL_TEST
1587 /* - - - - - - - - - */
1588
1589 /* Execute the linear XYZ device model */
domodel(double * out,double * in)1590 static void domodel(double *out, double *in) {
1591 double tmp[3];
1592 int i, j;
1593 double col[3][3]; /* sRGB additive colorant values in XYZ :- [out][in] */
1594
1595 col[0][0] = 0.412424; /* X from R */
1596 col[0][1] = 0.357579; /* X from G */
1597 col[0][2] = 0.180464; /* X from B */
1598 col[1][0] = 0.212656; /* Y from R */
1599 col[1][1] = 0.715158; /* Y from G */
1600 col[1][2] = 0.0721856; /* Y from B */
1601 col[2][0] = 0.0193324; /* Z from R */
1602 col[2][1] = 0.119193; /* Z from G */
1603 col[2][2] = 0.950444; /* Z from B */
1604
1605 #ifdef SPECIAL_TEST_GAMMA
1606 tmp[0] = pow(in[0], 1.9);
1607 tmp[1] = pow(in[1], 2.0);
1608 tmp[2] = pow(in[2], 2.1);
1609 #else
1610 tmp[0] = in[0];
1611 tmp[1] = in[1];
1612 tmp[2] = in[2];
1613 #endif
1614
1615 for (j = 0; j < 3; j++) {
1616 out[j] = 0.0;
1617 for (i = 0; i < 3; i++)
1618 out[j] += col[j][i] * tmp[i];
1619 }
1620 }
1621
1622 #ifdef SPECIAL_FORCE
1623 /* Function to pass to rspl to set nodes against */
1624 /* synthetic model. */
1625 static void
set_rspl_out1(void * pp,double * out,double * in)1626 set_rspl_out1(
1627 void *pp, /* relativectx structure */
1628 double *out, /* output value */
1629 double *in /* input value */
1630 ) {
1631 xfit *p = (xfit *)pp;
1632 double tt[3], tout[3];
1633
1634 /* Convert the input' values to input values */
1635 xfit_invinpscurves(p, tt, in);
1636
1637 /* Synthetic linear rgb->XYZ model */
1638 domodel(tout, tt);
1639
1640 /* Apply abs->rel white point adjustment */
1641 icmMulBy3x3(tout, p->cmat, tout);
1642
1643 #ifdef SPECIAL_TEST_LAB
1644 icmXYZ2Lab(&icmD50, tout, tout);
1645 #endif
1646 /* And then convert them back to clut values */
1647 xfit_invoutcurves(p, tout, tout);
1648
1649 #ifdef DEBUG
1650 printf("~1 changing %f %f %f -> %f %f %f\n", out[0], out[1], out[2], tout[0], tout[1], tout[2]);
1651 #endif
1652
1653 out[0] = tout[0];
1654 out[1] = tout[1];
1655 out[2] = tout[2];
1656 }
1657
1658 #endif /* SPECIAL_FORCE */
1659
1660 /* - - - - - - - - - */
1661 /* Do the fitting. */
1662 /* return nz on error */
1663 /* 1 = malloc or other error */
xfit_fit(struct _xfit * p,int flags,int di,int fdi,int rsplflags,double * wp,double * dw,double wpscale,double * dgw,cow * ipoints,int nodp,icxMatrixModel * skm,double in_min[MXDI],double in_max[MXDI],int gres[MXDI],double out_min[MXDO],double out_max[MXDO],double smooth,double oavgdev[MXDO],double demph,int iord[],int sord[],int oord[],double shp_smooth[MXDI],double out_smooth[MXDO],optcomb tcomb,void * cntx2,double (* to_de2)(void * cntx,double * in1,double * in2),double (* to_dde2)(void * cntx,double dout[2][MXDIDO],double * in1,double * in2))1664 static int xfit_fit(
1665 struct _xfit *p,
1666 int flags, /* Flag values */
1667 int di, /* Input dimensions */
1668 int fdi, /* Output dimensions */
1669 int rsplflags, /* clut rspl creation flags */
1670 double *wp, /* if flags & XFIT_OUT_WP_REL or XFIT_OUT_WP_REL_US, */
1671 /* Initial white point, returns final wp */
1672 double *dw, /* Device white value to adjust to be D50 */
1673 double wpscale, /* If >= 0.0 scale final wp */
1674 double *dgw, /* Device space gamut boundary white for XFIT_OUT_WP_REL_US */
1675 /* (ie. RGB 1,1,1 CMYK 0,0,0,0, etc) */
1676 cow *ipoints, /* Array of data points to fit - referece taken */
1677 int nodp, /* Number of data points */
1678 icxMatrixModel *skm, /* Optional skeleton model (used for input profiles) */
1679 double in_min[MXDI], /* Input value scaling/domain minimum */
1680 double in_max[MXDI], /* Input value scaling/domain maximum */
1681 int gres[MXDI], /* clut resolutions being optimised for/returned */
1682 double out_min[MXDO], /* Output value scaling/range minimum */
1683 double out_max[MXDO], /* Output value scaling/range maximum */
1684 // co *bpo, /* If != NULL, black point override in same spaces as ipoints */
1685 double smooth, /* clut rspl smoothing factor */
1686 double oavgdev[MXDO], /* Average output value deviation */
1687 double demph, /* dark emphasis factor for cLUT grid res. */
1688 int iord[], /* Order of input pos/shaper curve for each dimension */
1689 int sord[], /* Order of input sub-grid shaper curve (not used) */
1690 int oord[], /* Order of output shaper curve for each dimension */
1691 double shp_smooth[MXDI],/* Smoothing factors for each curve, nom = 1.0 */
1692 double out_smooth[MXDO],
1693 optcomb tcomb, /* Flag - target elements to fit. */
1694 void *cntx2, /* Context of callbacks */
1695 /* Callback to convert two fit values delta E squared */
1696 double (*to_de2)(void *cntx, double *in1, double *in2),
1697 /* Same as above, with partial derivatives */
1698 double (*to_dde2)(void *cntx, double dout[2][MXDIDO], double *in1, double *in2)
1699 ) {
1700 int i, e, f;
1701 double *b; /* Base of parameters for this section */
1702 int poff;
1703
1704 double powtol = POWTOL1; /* powell/conjgrad initial tollerance */
1705 int maxits = MAXITS1; /* powell/conjgrad initial maximum itterations */
1706
1707 if (tcomb & oc_io) /* If we're doing anything, we need the matrix */
1708 tcomb |= oc_m;
1709
1710 p->flags = flags;
1711 if (flags & XFIT_VERB)
1712 p->verb = 1;
1713 else
1714 p->verb = 0;
1715 p->di = di;
1716 p->fdi = fdi;
1717 p->wp = wp; /* Take reference, so modified wp can be returned */
1718 p->dw = dw;
1719 p->nodp = nodp;
1720 p->skm = skm; /* This isn't current used by profin, because it doesn't help.. */
1721 p->ipoints = ipoints;
1722 p->tcomb = tcomb;
1723 p->cntx2 = cntx2;
1724 for (e = 0; e < di; e++)
1725 p->gres[e] = gres[e];
1726 p->to_de2 = to_de2;
1727 p->to_dde2 = to_dde2;
1728
1729 #ifdef DEBUG
1730 printf("xfit_fit called with flags = 0x%x, di = %d, fdi = %d, nodp = %d, tcomb = 0x%x\n",flags,di,fdi,nodp,tcomb);
1731 #endif
1732
1733 //printf("~1 out min = %f %f %f max = %f %f %f\n", out_min[0], out_min[1], out_min[2], out_max[0], out_max[1], out_max[2]);
1734
1735 /* Sanity protect shaper orders */
1736 /* and save scaling and smoothness factors. */
1737 p->sm_iluord = MXLUORD+1;
1738 for (e = 0; e < di; e++) {
1739 if (iord[e] > MXLUORD)
1740 p->iluord[e] = MXLUORD;
1741 else
1742 p->iluord[e] = iord[e];
1743 if (p->iluord[e] < p->sm_iluord)
1744 p->sm_iluord = p->iluord[e];
1745 p->in_min[e] = in_min[e];
1746 p->in_max[e] = in_max[e];
1747 p->shp_smooth[e] = shp_smooth[e];
1748 }
1749 for (f = 0; f < fdi; f++) {
1750 if (oord[f] > MXLUORD)
1751 p->oluord[f] = MXLUORD;
1752 else
1753 p->oluord[f] = oord[f];
1754 p->out_min[f] = out_min[f];
1755 p->out_max[f] = out_max[f];
1756 p->out_smooth[f] = out_smooth[f];
1757 }
1758
1759
1760 /* Compute parameter offset and count information */
1761 p->shp_off = 0;
1762 for (poff = p->shp_off, p->shp_cnt = 0, e = 0; e < di; e++) {
1763 p->shp_offs[e] = poff;
1764 p->shp_cnt += p->iluord[e];
1765 poff += p->iluord[e];
1766 }
1767
1768 p->mat_off = p->shp_off + p->shp_cnt;
1769 for (poff = p->mat_off, p->mat_cnt = 0, f = 0; f < fdi; f++) {
1770 p->mat_offs[f] = poff;
1771 p->mat_cnt += (1 << di);
1772 poff += (1 << di);
1773 }
1774
1775 p->out_off = p->mat_off + p->mat_cnt;
1776 for (poff = p->out_off, p->out_cnt = 0, f = 0; f < fdi; f++) {
1777 p->out_offs[f] = poff;
1778 p->out_cnt += p->oluord[f];
1779 poff += p->oluord[f];
1780 }
1781
1782 p->pos_off = p->out_off + p->out_cnt;
1783 for (poff = p->pos_off, p->pos_cnt = 0, e = 0; e < di; e++) {
1784 p->pos_offs[e] = poff;
1785 p->pos_cnt += p->iluord[e];
1786 poff += p->iluord[e];
1787 }
1788
1789 p->tot_cnt = p->shp_cnt + p->mat_cnt + p->out_cnt + p->pos_cnt;
1790 if (p->tot_cnt > MXPARMS)
1791 error("xfit_fit: assert tot_cnt exceeds MXPARMS");
1792
1793 /* Allocate space for parameter values */
1794 if (p->v != NULL) {
1795 free(p->v);
1796 p->v = NULL;
1797 }
1798 if (p->wv != NULL) {
1799 free(p->wv);
1800 p->wv = NULL;
1801 }
1802 if (p->sa != NULL) {
1803 free(p->sa);
1804 p->sa = NULL;
1805 }
1806 if ((p->v = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
1807 return 1;
1808 if ((p->wv = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
1809 return 1;
1810 if ((p->sa = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
1811 return 1;
1812
1813 /* Setup initial white point abs->rel conversions */
1814 if ((p->flags & XFIT_OUT_WP_REL) != 0) {
1815 icmXYZNumber _wp;
1816
1817 icmAry2XYZ(_wp, p->wp);
1818
1819 /* Absolute->Aprox. Relative Adaptation matrix and */
1820 /* Aproximate relative to absolute conversion matrix */
1821 if (p->picc != NULL) {
1822 p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, p->toAbs, p->fromAbs, icmD50, _wp);
1823 } else {
1824 icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
1825 icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
1826 }
1827
1828 if (p->verb) {
1829 double lab[3];
1830 icmXYZ2Lab(&icmD50, lab, p->wp);
1831 printf("Initial White Point XYZ %f %f %f, Lab %f %f %f\n",
1832 p->wp[0], p->wp[1], p->wp[2], lab[0], lab[1], lab[2]);
1833 }
1834 } else {
1835 icmSetUnity3x3(p->fromAbs);
1836 icmSetUnity3x3(p->toAbs);
1837 }
1838
1839 /* Setup input position/shape curves to be linear initially */
1840 b = p->v + p->shp_off;
1841 for (e = 0; e < di; b += p->iluord[e], e++) {
1842 for (i = 0; i < p->iluord[e]; i++) {
1843 b[i] = 0.0;
1844 }
1845 }
1846
1847 /* Setup matrix to be pure colorant' values initially */
1848 b = p->v + p->mat_off;
1849 for (e = 0; e < (1 << di); e++) { /* For each colorant combination */
1850 int j, k, bk = 0;
1851 double bdif = 1e6;
1852 double ov[MXDO];
1853
1854 /* Search the patch list to find the one closest to this input combination */
1855 for (k = 0; k < p->nodp; k++) {
1856 double dif = 0.0;
1857 for (j = 0; j < di; j++) {
1858 double tt;
1859 if (e & (1 << j))
1860 tt = p->in_max[j] - p->ipoints[k].p[j];
1861 else
1862 tt = p->in_min[j] - p->ipoints[k].p[j];
1863 dif += tt * tt;
1864 }
1865 if (dif < bdif) { /* best so far */
1866 bdif = dif;
1867 bk = k;
1868 if (dif < 0.001)
1869 break; /* Don't bother looking further */
1870 }
1871 }
1872 xfit_abs_to_rel(p, ov, p->ipoints[bk].v);
1873
1874 for (f = 0; f < fdi; f++)
1875 b[f * (1 << di) + e] = ov[f];
1876 }
1877
1878 /* Setup output curves to be linear initially */
1879 b = p->v + p->out_off;
1880 for (f = 0; f < fdi; b += p->oluord[f], f++) {
1881 for (i = 0; i < p->oluord[f]; i++) {
1882 b[i] = 0.0;
1883 }
1884 }
1885
1886 /* Setup positioning curves to be linear initially */
1887 b = p->v + p->pos_off;
1888 for (e = 0; e < di; b += p->iluord[e], e++) {
1889 for (i = 0; i < p->iluord[e]; i++) {
1890 b[i] = 0.0;
1891 }
1892 }
1893
1894 /* Create copy of input points with output converted to white relative */
1895 if (p->rpoints == NULL) {
1896 if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
1897 return 1;
1898 }
1899 for (i = 0; i < p->nodp; i++) {
1900 p->rpoints[i].w = p->ipoints[i].w;
1901 for (e = 0; e < di; e++)
1902 p->rpoints[i].p[e] = p->ipoints[i].p[e];
1903 for (f = 0; f < fdi; f++)
1904 p->rpoints[i].v[f] = p->ipoints[i].v[f];
1905
1906 /* out -> rout */
1907 xfit_abs_to_rel(p, p->rpoints[i].v, p->rpoints[i].v);
1908 }
1909
1910 /* Allocate array of pseudo-inverse matricies */
1911 if ((p->flags & XFIT_FM_INPUT) != 0 && p->piv == NULL) {
1912 if ((p->piv = (xfit_piv *)malloc(p->nodp * sizeof(xfit_piv))) == NULL)
1913 return 1;
1914 }
1915
1916 /* Allocate array of span DE's for current opt channel */
1917 {
1918 int lres = 0;
1919 for (e = 0; e < di; e++) {
1920 if (p->gres[e] > lres)
1921 lres = p->gres[e];
1922 }
1923 if ((p->uerrv = (double *)malloc(lres * sizeof(double))) == NULL)
1924 return 1;
1925 }
1926
1927 /* Do the fitting one part at a time, then together */
1928 /* Shaper curves are created if position or shaper curves are requested */
1929
1930 /* Fit just the matrix */
1931 if ((p->tcomb & oc_ipo) != 0
1932 && (p->tcomb & oc_m) == oc_m) { /* Only bother with matrix if in and/or out */
1933 double rerr;
1934
1935 if (p->verb)
1936 printf("About to optimise temporary matrix\n");
1937
1938 /* Setup pseudo-inverse if we need it */
1939 if (p->flags & XFIT_FM_INPUT)
1940 setup_piv(p);
1941
1942 #ifdef DEBUG
1943 printf("\nBefore matrix opt:\n");
1944 dump_xfit(p);
1945 #endif
1946 /* Optimise matrix on its own */
1947 p->opt_ssch = 0;
1948 p->opt_ch = -1;
1949 p->opt_msk = oc_m;
1950 setup_xfit(p, p->wv, p->sa, 0.0, 0.5);
1951
1952 #ifdef NODDV
1953 if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
1954 xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
1955 warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
1956 #else
1957 if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
1958 xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
1959 warning("xfit_fit: Conjgrad failed to converge, residual error = %f", rerr);
1960 #endif
1961 for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
1962 p->v[p->opt_off + i] = p->wv[i];
1963
1964 #ifdef DEBUG
1965 printf("\nAfter matrix opt:\n");
1966 dump_xfit(p);
1967
1968 #endif
1969 }
1970
1971 /* Optimise input and matrix together */
1972 if ((p->tcomb & oc_im) == oc_im) {
1973 double rerr;
1974 int sm_iluord = p->sm_iluord;
1975
1976 if (p->verb)
1977 printf("About to optimise a common ord 0 input curve and matrix\n");
1978
1979 /* Setup pseudo-inverse if we need it */
1980 if (p->flags & XFIT_FM_INPUT)
1981 setup_piv(p);
1982
1983 p->opt_ssch = 1;
1984 p->sm_iluord = 1; /* Do a single order for first up */
1985 p->opt_ch = -1;
1986 p->opt_msk = oc_im;
1987 setup_xfit(p, p->wv, p->sa, 0.5, 0.3);
1988
1989 #ifdef NODDV
1990 if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
1991 xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
1992 #ifdef DEBUG
1993 warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
1994 #endif
1995 }
1996 #else
1997 if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
1998 xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
1999 #ifdef DEBUG
2000 warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
2001 #endif
2002 }
2003 #endif /* !NODDV */
2004 for (e = 0; e < di; e++) { /* Copy optimised values back */
2005 for (i = 0; i < p->sm_iluord; i++)
2006 p->v[p->shp_offs[e] + i] = p->wv[i];
2007 for (; i < p->iluord[e]; i++)
2008 p->v[p->shp_offs[e] + i] = 0.0;
2009 }
2010 for (i = p->sm_iluord; i < p->opt_cnt; i++)
2011 p->v[p->mat_off + i - p->sm_iluord] = p->wv[i];
2012 #ifdef DEBUG
2013 printf("\nAfter single input and matrix opt:\n");
2014 dump_xfit(p);
2015 #endif
2016
2017 /* - - - - - - - - - - - */
2018 if (p->verb)
2019 printf("About to optimise a common input curve and matrix\n");
2020
2021 /* Setup pseudo-inverse if we need it */
2022 if (p->flags & XFIT_FM_INPUT)
2023 setup_piv(p);
2024
2025 p->opt_ssch = 1;
2026 p->sm_iluord = sm_iluord; /* restore this */
2027 p->opt_ch = -1;
2028 p->opt_msk = oc_im;
2029 setup_xfit(p, p->wv, p->sa, 0.5, 0.3);
2030
2031 #ifdef NODDV
2032 if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
2033 xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
2034 #ifdef DEBUG
2035 warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
2036 #endif
2037 }
2038 #else
2039 if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
2040 xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
2041 #ifdef DEBUG
2042 warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
2043 #endif
2044 }
2045 #endif /* !NODDV */
2046 for (e = 0; e < di; e++) { /* Copy optimised values back */
2047 for (i = 0; i < p->sm_iluord; i++)
2048 p->v[p->shp_offs[e] + i] = p->wv[i];
2049 for (; i < p->iluord[e]; i++)
2050 p->v[p->shp_offs[e] + i] = 0.0;
2051 }
2052 for (i = p->sm_iluord; i < p->opt_cnt; i++)
2053 p->v[p->mat_off + i - p->sm_iluord] = p->wv[i];
2054 #ifdef DEBUG
2055 printf("\nAfter single input and matrix opt:\n");
2056 dump_xfit(p);
2057 #endif
2058
2059 /* - - - - - - - - - - - */
2060 if (p->verb)
2061 printf("About to optimise input curves and matrix\n");
2062
2063 if ((p->tcomb & oc_mo) != oc_mo) { /* If this will be last fit */
2064 powtol = POWTOL;
2065 maxits = MAXITS;
2066 }
2067
2068 /* Setup pseudo-inverse if we need it */
2069 if (p->flags & XFIT_FM_INPUT)
2070 setup_piv(p);
2071
2072 p->opt_ssch = 0;
2073 p->opt_ch = -1;
2074 p->opt_msk = oc_im;
2075 setup_xfit(p, p->wv, p->sa, 0.5, 0.3);
2076 /* Suppress the warnings the first time through - it's better to cut off the */
2077 /* itterations and move on to the output curve, and worry about it not */
2078 /* converging the second time through. */
2079 #ifdef NODDV
2080 if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
2081 xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
2082 #ifdef DEBUG
2083 warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
2084 #endif
2085 }
2086 #else
2087 if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
2088 xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
2089 #ifdef DEBUG
2090 warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
2091 #endif
2092 }
2093 #endif /* !NODDV */
2094 for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
2095 p->v[p->opt_off + i] = p->wv[i];
2096 #ifdef DEBUG
2097 printf("\nAfter input and matrix opt:\n");
2098 dump_xfit(p);
2099 #endif
2100 }
2101
2102 /* Optimise the matrix and output curves together */
2103 if ((p->tcomb & oc_mo) == oc_mo) {
2104 double rerr;
2105
2106 if (p->verb)
2107 printf("About to optimise output curves and matrix\n");
2108
2109 if ((p->tcomb & oc_im) != oc_im) { /* If this will be last fit */
2110 powtol = POWTOL;
2111 maxits = MAXITS;
2112 }
2113
2114 /* Setup pseudo-inverse if we need it */
2115 if (p->flags & XFIT_FM_INPUT)
2116 setup_piv(p);
2117
2118 p->opt_ssch = 0;
2119 p->opt_ch = -1;
2120 p->opt_msk = oc_mo;
2121 setup_xfit(p, p->wv, p->sa, 0.3, 0.3);
2122 #ifdef NODDV
2123 if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
2124 xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
2125 warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
2126 #else
2127 if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits, xfitfunc,
2128 dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
2129 warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
2130 #endif
2131 for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
2132 p->v[p->opt_off + i] = p->wv[i];
2133 #ifdef DEBUG
2134 printf("\nAfter output opt:\n");
2135 dump_xfit(p);
2136 #endif
2137
2138 /* Optimise input and matrix together again, after altering matrix */
2139 if ((p->tcomb & oc_im) == oc_im) {
2140
2141 if (p->verb)
2142 printf("About to optimise input curves and matrix again\n");
2143
2144
2145 #ifndef NODDV
2146 if ((p->tcomb & oc_imo) != oc_imo) /* If this will be last fit */
2147 #endif
2148 {
2149 powtol = POWTOL;
2150 maxits = MAXITS;
2151 }
2152
2153 /* Setup pseudo-inverse if we need it */
2154 if (p->flags & XFIT_FM_INPUT)
2155 setup_piv(p);
2156
2157 p->opt_ssch = 0;
2158 p->opt_ch = -1;
2159 p->opt_msk = oc_im;
2160 setup_xfit(p, p->wv, p->sa, 0.2, 0.2);
2161 #ifdef NODDV
2162 if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
2163 xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
2164 warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
2165 #else
2166 if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
2167 xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
2168 warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
2169 #endif
2170 for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
2171 p->v[p->opt_off + i] = p->wv[i];
2172 }
2173 #ifdef DEBUG
2174 printf("\nAfter 2nd input and matrix opt:\n");
2175 dump_xfit(p);
2176 #endif
2177
2178 #ifndef NODDV
2179 /* Optimise all together */
2180 /* (This is very slow using powell) */
2181 if ((p->tcomb & oc_imo) == oc_imo) {
2182
2183 if (p->verb)
2184 printf("About to optimise input, matrix and output together\n");
2185
2186 /* Setup pseudo-inverse if we need it */
2187 if (p->flags & XFIT_FM_INPUT)
2188 setup_piv(p);
2189
2190 p->opt_ssch = 0;
2191 p->opt_ch = -1;
2192 p->opt_msk = oc_imo;
2193 setup_xfit(p, p->wv, p->sa, 0.1, 0.1);
2194 if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
2195 xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
2196 warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
2197 for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
2198 p->v[p->opt_off + i] = p->wv[i];
2199
2200 /* Setup final pseudo-inverse from shaper/matrix/out */
2201 if (p->flags & XFIT_FM_INPUT)
2202 setup_piv(p);
2203 }
2204
2205 #ifdef DEBUG
2206 printf("\nAfter all together opt:\n");
2207 dump_xfit(p);
2208 #endif
2209
2210 #endif /* !NODDV */
2211
2212 /* Adjust output curve white point. */
2213 /* This is for the benefit of the B2A table */
2214 if (p->flags & XFIT_OUT_ZERO) {
2215
2216 if (p->verb)
2217 printf("About to adjust a and b output curves for white point\n");
2218
2219 for (f = 1; f < 3 && f < p->fdi; f++) {
2220 p->opt_ch = f;
2221 p->wv[0] = p->v[p->out_offs[f]]; /* Current parameter value */
2222 p->sa[0] = 0.1; /* Search radius */
2223 if (powell(&rerr, 1, p->wv, p->sa, 0.0000001, 1000,
2224 symoptfunc, (void *)p, NULL, NULL) != 0)
2225 error("xfit_fit: Powell failed to converge, residual error = %f",rerr);
2226 p->v[p->out_offs[f]] = p->wv[0]; /* Copy results back */
2227 }
2228 }
2229 }
2230
2231 /* In case we don't generate position curves, */
2232 /* copy the input curves to the position, so that */
2233 /* ipos is computed correctly */
2234 {
2235 double *bb;
2236
2237 b = p->v + p->shp_off;
2238 bb = p->v + p->pos_off;
2239 for (e = 0; e < di; b += p->iluord[e], bb += p->iluord[e], e++) {
2240 for (i = 0; i < p->iluord[e]; i++) {
2241 bb[i] = b[i];
2242 }
2243 }
2244 }
2245
2246
2247 /* If we want position curves, generate them */
2248 /* (This could possibly be improved by using some sort */
2249 /* of optimization drivel approach rather than the predictive */
2250 /* method used here.) */
2251 if (p->tcomb & oc_p) {
2252 int ee;
2253
2254 if (p->verb)
2255 printf("About to create grid position input curves\n");
2256
2257 /* Allocate in->rout duplicate point set */
2258 if (p->rpoints == NULL) {
2259 if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
2260 return 1;
2261 }
2262
2263 /* Create a set of 1D rspl setup points that contains */
2264 /* the residual error */
2265 for (i = 0; i < p->nodp; i++) {
2266 double tv[MXDO]; /* Target output value */
2267 double mv[MXDO]; /* Model output value */
2268 double ev;
2269
2270 xfit_abs_to_rel(p, tv, p->ipoints[i].v);
2271 xfit_shmatsh(p, mv, p->ipoints[i].p);
2272
2273 /* Evaluate the error squared */
2274 if (p->flags & XFIT_FM_INPUT) {
2275 double pp[MXDI];
2276 for (e = 0; e < di; e++)
2277 pp[e] = p->ipoints[i].p[e];
2278 for (f = 0; f < fdi; f++) {
2279 double t1 = tv[f] - mv[f]; /* Error in output */
2280
2281 /* Create input point offset by equivalent delta to output point */
2282 /* error, in proportion to the partial derivatives for that output. */
2283 for (e = 0; e < di; e++)
2284 pp[e] += t1 * p->piv[i].ide[f][e];
2285 }
2286 ev = p->to_de2(p->cntx2, pp, p->ipoints[i].p);
2287 } else {
2288 ev = p->to_de2(p->cntx2, mv, tv);
2289 }
2290 p->rpoints[i].v[0] = ev;
2291 p->rpoints[i].w = p->ipoints[i].w;
2292 }
2293
2294 /* Do each input axis in turn */
2295 for (ee = 0; ee < p->di; ee++) {
2296 rspl *resid;
2297 double imin[1],imax[1],omin[1],omax[1];
2298 int resres[1] = { 1024 };
2299 #define NPGP 100
2300 mcv *posc;
2301 mcvco pgp[NPGP];
2302 double vo, vs;
2303 double *pms;
2304
2305 /* Create a rspl that plots the residual error */
2306 /* vs the axis value */
2307 for (i = 0; i < p->nodp; i++)
2308 p->rpoints[i].p[0] = p->ipoints[i].p[ee];
2309
2310 imin[0] = in_min[ee];
2311 imax[0] = in_max[ee];
2312 omin[0] = 0.0;
2313 omax[0] = 0.0;
2314
2315 if ((resid = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL)
2316 return 1;
2317
2318 resid->fit_rspl_w(resid, RSPLFLAGS, p->rpoints, p->nodp, imin, imax, resres,
2319 omin, omax, 2.0, NULL, NULL);
2320
2321 #ifdef DEBUG_PLOT
2322 {
2323 #define XRES 100
2324 double xx[XRES];
2325 double y1[XRES];
2326
2327 printf("Input residual error channel %d\n",ee);
2328 for (i = 0; i < XRES; i++) {
2329 co pp;
2330 double x;
2331 x = i/(double)(XRES-1);
2332 xx[i] = x = x * (imax[0] - imin[0]) + imin[0];
2333 pp.p[0] = xx[i];
2334 resid->interp(resid, &pp);
2335 y1[i] = pp.v[0];
2336 if (y1[i] < 0.0)
2337 y1[i] = 0.0;
2338 y1[i] = pow(y1[i], 0.5); /* Convert from error^2 to error */
2339 }
2340 do_plot(xx,y1,NULL,NULL,XRES);
2341 }
2342 #endif /* DEBUG_PLOT */
2343
2344 /* Create a set of guide points that contain */
2345 /* the accumulated residual error vs. the axis value */
2346 for (i = 0; i < NPGP; i++) {
2347 co pp;
2348 double vv;
2349
2350 pp.p[0] = i/(NPGP-1.0);
2351 resid->interp(resid, &pp);
2352
2353 pgp[i].p = (pp.p[0] - in_min[ee])/(in_max[ee] - in_min[ee]);
2354 pgp[i].w = 1.0;
2355
2356 vv = pp.v[0];
2357 if (vv < 0.0)
2358 vv = 0.0;
2359 vv = pow(vv, 0.5); /* Convert from error^2 to error */
2360 vv += PSHAPE_MINE; /* In case error is near zero */
2361 vv = pow(vv, PSHAPE_DIST); /* Agressivness of grid distribution */
2362
2363 if (i == 0)
2364 pgp[i].v = vv;
2365 else
2366 pgp[i].v = pgp[i-1].v + vv;
2367
2368 }
2369 resid->del(resid);
2370
2371 /* Normalize the output range */
2372 vo = pgp[0].v;
2373 vs = pgp[NPGP-1].v - vo;
2374
2375 for (i = 0; i < NPGP; i++) {
2376 pgp[i].v = (pgp[i].v - vo)/vs;
2377
2378 /* Apply any dark emphasis */
2379 if (demph > 1.0) {
2380 pgp[i].v = icx_powlike(pgp[i].v, 1.0/demph);
2381 }
2382 }
2383 /* Fit the non-monotonic parameters to the guide points */
2384 if ((posc = new_mcv_noos()) == NULL)
2385 return 1;
2386
2387 posc->fit(posc, 0, p->iluord[ee], pgp, NPGP, 0.1); // ~~99
2388
2389 #ifdef DEBUG_PLOT
2390 {
2391 #define XRES 100
2392 double xx[XRES];
2393 double y1[XRES];
2394
2395 printf("Position curve %d\n",ee);
2396 for (i = 0; i < XRES; i++) {
2397 xx[i] = i/(double)(XRES-1);
2398 y1[i] = posc->interp(posc, xx[i]);
2399 }
2400 do_plot(xx,y1,NULL,NULL,XRES);
2401 }
2402 #endif /* DEBUG_PLOT */
2403
2404 /* Transfer parameters to xfit pos (skip offset and scale) */
2405 posc->get_params(posc, &pms);
2406
2407 for (i = 0; i < p->iluord[ee]; i++) {
2408 p->v[p->pos_offs[ee] + i] = pms[i+2];
2409 }
2410 //p->v[p->in_offs[ee]] = -1.5;
2411
2412 free(pms);
2413 posc->del(posc);
2414 }
2415 }
2416
2417 #ifdef DEBUG
2418 printf("Final parameters:\n");
2419 dump_xfit(p);
2420 #endif
2421
2422 #ifdef DEBUG_PLOT
2423 {
2424 #define XRES 100
2425 double xx[XRES];
2426 double y1[XRES];
2427
2428 for (e = 0; e < p->di; e++) {
2429 printf("Input position curve channel %d\n",e);
2430 for (i = 0; i < XRES; i++) {
2431 double x;
2432 x = i/(double)(XRES-1);
2433 xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
2434 y1[i] = xfit_poscurve(p, x, e);
2435 }
2436 do_plot(xx,y1,NULL,NULL,XRES);
2437 }
2438
2439 for (e = 0; e < p->di; e++) {
2440 printf("Input shape curve channel %d\n",e);
2441 for (i = 0; i < XRES; i++) {
2442 double x;
2443 x = i/(double)(XRES-1);
2444 xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
2445 y1[i] = xfit_shpcurve(p, x, e);
2446 }
2447 do_plot(xx,y1,NULL,NULL,XRES);
2448 }
2449
2450 for (e = 0; e < p->di; e++) {
2451 printf("Combined input curve channel %d\n",e);
2452 for (i = 0; i < XRES; i++) {
2453 double x;
2454 x = i/(double)(XRES-1);
2455 xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
2456 y1[i] = p->incurve(p, x, e);
2457 }
2458 do_plot(xx,y1,NULL,NULL,XRES);
2459 }
2460
2461 for (f = 0; f < p->fdi; f++) {
2462 printf("Output curve channel %d\n",f);
2463 for (i = 0; i < XRES; i++) {
2464 double x;
2465 x = i/(double)(XRES-1);
2466 xx[i] = x = x * (p->out_max[f] - p->out_min[f]) + p->out_min[f];
2467 y1[i] = p->outcurve(p, x, f);
2468 }
2469 do_plot(xx,y1,NULL,NULL,XRES);
2470 }
2471 }
2472 #endif /* DEBUG_PLOT */
2473
2474 /* Create final clut rspl using the established pos/shape/output curves */
2475 /* and white point */
2476 if (flags & XFIT_MAKE_CLUT) {
2477 double *ipos[MXDI];
2478
2479 /* Create an in' -> rout' scattered test point set */
2480 if (p->rpoints == NULL) {
2481 if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
2482 return 1;
2483 }
2484 for (i = 0; i < p->nodp; i++) {
2485 p->rpoints[i].w = p->ipoints[i].w;
2486 xfit_inpscurves(p, p->rpoints[i].p, p->ipoints[i].p);
2487 for (f = 0; f < fdi; f++)
2488 p->rpoints[i].v[f] = p->ipoints[i].v[f];
2489 xfit_abs_to_rel(p, p->rpoints[i].v, p->rpoints[i].v);
2490 xfit_invoutcurves(p, p->rpoints[i].v, p->rpoints[i].v);
2491
2492 if (p->skm) {
2493 /* Look up the skeleton value */
2494 double skval[MXDO];
2495 p->skm->lookup(p->skm, skval, p->ipoints[i].p);
2496 xfit_abs_to_rel(p, skval, skval);
2497 xfit_invoutcurves(p, skval, skval);
2498 //printf("~1 point %d at %f %f %f, targ %f %f %f skm %f %f %f\n", i,p->ipoints[i].p[0],p->ipoints[i].p[1],p->ipoints[i].p[2], p->rpoints[i].v[0],p->rpoints[i].v[1],p->rpoints[i].v[2], skval[0], skval[1], skval[2]);
2499 /* Subtract it from value at this point, */
2500 /* so rspl will fit difference to skeleton model */
2501 for (f = 0; f < fdi; f++)
2502 p->rpoints[i].v[f] -= skval[f];
2503 }
2504 //printf("~1 point %d, w %f, %f %f %f %f -> %f %f %f\n",i,p->rpoints[i].w,p->rpoints[i].p[0], p->rpoints[i].p[1], p->rpoints[i].p[2], p->rpoints[i].p[3],p->rpoints[i].v[0], p->rpoints[i].v[1], p->rpoints[i].v[2]);
2505 }
2506
2507 /* Create ipos[] arrays, that hold the shaper space */
2508 /* grid position due to the positioning curves. */
2509 /* This tells the rspl scattered data interpolator */
2510 /* the grid spacing that smoothness should be */
2511 /* measured against. */
2512 #ifdef DEBUG
2513 printf("~1 about to setup ipos\n");
2514 #endif
2515 for (e = 0; e < p->di; e++) {
2516 //printf("~1 e = %d\n",e);
2517 if ((ipos[e] = (double *)malloc((p->gres[e]) * sizeof(double))) == NULL)
2518 return 1;
2519 //printf("~1 about to do %d spans\n",p->gres[e]);
2520 for (i = 0; i < p->gres[e]; i++) {
2521 double cv;
2522 cv = (double)i/p->gres[e];
2523
2524 //printf("~1 i = %d, pos space = %f\n",i,cv);
2525 /* Inverse lookup grid position through positioning curve */
2526 /* to give device space type value */
2527 cv = icxInvTransFunc(p->v + p->pos_offs[e], p->iluord[e], cv);
2528
2529 //printf("~1 dev space = %f\n",cv);
2530 /* Forward lookup device type value through the shaper curve */
2531 /* to give value in shape linearized space. */
2532 cv = icxTransFunc(p->v + p->shp_offs[e], p->iluord[e], cv);
2533 //printf("~1 shape space = %f\n",cv);
2534 ipos[e][i] = cv;
2535 #ifdef DEBUG
2536 printf("~1 ipos[%d][%d] = %f\n",e,i,cv);
2537 #endif
2538 }
2539 }
2540
2541 if (p->clut != NULL)
2542 p->clut->del(p->clut);
2543 if ((p->clut = new_rspl(RSPL_NOFLAGS, di, fdi)) == NULL)
2544 return 1;
2545
2546 if (p->verb)
2547 printf("Create final clut from scattered data\n");
2548
2549 #ifdef EXTEND_GRID
2550 #define XN EXTEND_GRID_BYN
2551 /* Try increasing the grid by one row all around */
2552 {
2553 #pragma message("!!!!!!!!!!!! Experimental rspl fitting resolution !!!!!!!!!")
2554 double xin_min[MXDI];
2555 double xin_max[MXDI];
2556 int xgres[MXDI];
2557 double del;
2558 double *xipos[MXDI];
2559
2560 for (e = 0; e < p->di; e++) {
2561 del = (in_max[e] - in_min[e])/(gres[e]-1.0); /* Extension */
2562 xin_min[e] = in_min[e] - XN * del;
2563 xin_max[e] = in_max[e] + XN * del;
2564 xgres[e] = gres[e] + 2 * XN;
2565 //printf("~1 xgres %d, gres %d\n",xgres[e], gres[e]);
2566
2567 if ((xipos[e] = (double *)malloc((xgres[e]) * sizeof(double))) == NULL)
2568 return 1;
2569 for (i = 0; i < xgres[e]; i++) {
2570 if (i < XN) { /* Extrapolate bottom */
2571 xipos[e][i] = ipos[e][0] - (XN - i) * (ipos[e][1] - ipos[e][0]);
2572 //printf("~1 xipos[%d] %f from ipos[%d] %f and ipos[%d] %f\n",i,xipos[e][i],0,ipos[e][0],1,ipos[e][1]);
2573 } else if (i >= (xgres[e]-XN)) { /* Extrapolate top */
2574 xipos[e][i] = ipos[e][gres[e]-1] + (i - xgres[e] + XN + 1) * (ipos[e][gres[e]-1] - ipos[e][gres[e]-2]);
2575 //printf("~1 xipos[%d] %f from ipos[%d] %f and ipos[%d] %f\n",i,xipos[e][i],gres[e]-1,ipos[e][gres[e]-1],gres[e]-2,ipos[e][gres[e]-2]);
2576 } else {
2577 xipos[e][i] = ipos[e][i-XN];
2578 //printf("~1 xipos[%d] %f from ipos[%d] %f\n",i,xipos[e][i],i-XN,ipos[e][i-XN]);
2579 }
2580 }
2581 }
2582
2583 p->clut->fit_rspl_w(p->clut, rsplflags, p->rpoints, p->nodp, xin_min, xin_max, xgres,
2584 out_min, out_max, smooth, oavgdev, xipos);
2585
2586 for (e = 0; e < p->di; e++) {
2587 free(xipos[e]);
2588 }
2589 }
2590 #undef XN
2591 #else
2592 // if (p->skm) {
2593 // /* This doesn't seem to work as well as some explicit neutral axis points.. */
2594 // p->clut->fit_rspl_w_df(p->clut, rsplflags, p->rpoints, p->nodp, in_min, in_max, gres,
2595 // out_min, out_max, smooth, oavgdev, ipos, 1.0, (void *)p, skm_weak);
2596 // } else
2597 /* Normal multi-d scattered point fitting */
2598 p->clut->fit_rspl_w(p->clut, rsplflags, p->rpoints, p->nodp, in_min, in_max, gres,
2599 out_min, out_max, smooth, oavgdev, ipos);
2600 #endif
2601 if (p->verb)
2602 printf("\n");
2603
2604 for (e = 0; e < p->di; e++)
2605 free(ipos[e]);
2606
2607 /* If we used a skeleton model, add it into the resulting rspl values */
2608 if (p->skm) {
2609
2610 /* Undo the input point change to allow diagnostic code to work */
2611 for (i = 0; i < p->nodp; i++) {
2612 /* Look up the skeleton value */
2613 double skval[MXDO];
2614 p->skm->lookup(p->skm, skval, p->ipoints[i].p);
2615 xfit_abs_to_rel(p, skval, skval);
2616 xfit_invoutcurves(p, skval, skval);
2617
2618 /* Subtract it from value at this point, */
2619 /* so rspl will fit difference to skeleton model */
2620 for (f = 0; f < fdi; f++)
2621 p->rpoints[i].v[f] += skval[f];
2622 }
2623
2624 /* Undo the skm from the resultant rspl */
2625 p->clut->re_set_rspl(
2626 p->clut, /* this */
2627 0, /* Combination of flags */
2628 (void *)p, /* Opaque function context */
2629 skm_rspl_out /* Function to set from */
2630 );
2631
2632 }
2633
2634 /* The overall device to absolute conversion is now what we want */
2635 /* (as dictated by the points, weighting and best fit), */
2636 /* but we need to adjust the device to relative conversion */
2637 /* to make device white map exactly to D50, without touching */
2638 /* the overall absolute behaviour. */
2639 if (p->flags & XFIT_OUT_WP_REL) {
2640 co wcc; /* device white + aprox rel. white */
2641 icmXYZNumber _wp; /* Uncorrected dw maps to _wp */
2642
2643 if (p->verb)
2644 printf("Doing White point fine tune:\n");
2645
2646 /* See what the relative and absolute white point has turned out to be, */
2647 /* by looking up the device white in the current conversion */
2648 xfit_inpscurves(p, wcc.p, dw);
2649 p->clut->interp(p->clut, &wcc);
2650 xfit_outcurves(p, wcc.v, wcc.v);
2651 if (p->flags & XFIT_OUT_LAB)
2652 icmLab2XYZ(&icmD50, wcc.v, wcc.v);
2653
2654 if (p->verb) {
2655 double labwp[3];
2656 icmXYZ2Lab(&icmD50, labwp, wcc.v);
2657 printf("Before fine tune, rel WP = XYZ %f %f %f, Lab %f %f %f\n",
2658 wcc.v[0], wcc.v[1],wcc.v[2], labwp[0], labwp[1], labwp[2]);
2659 }
2660
2661 /* Matrix needed to correct approx rel wp to target D50 */
2662 icmAry2XYZ(_wp, wcc.v); /* Aprox relative target white point */
2663 if (p->picc != NULL) /* Correction */
2664 p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, NULL, p->cmat, icmD50, _wp);
2665 else
2666 icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->cmat);
2667
2668 /* Compute the actual white point, and return it to caller */
2669 icmMulBy3x3(wp, p->toAbs, wcc.v);
2670 icmAry2Ary(p->wp, wp);
2671
2672 /* Apply correction to fine tune rspl data. */
2673 /* NOTE: this doesn't always give us a perfect D50 white for */
2674 /* Lab PCS input profiles because the dev white may land */
2675 /* within a cell, and the clipping of Lab PCS values in the grid */
2676 /* may introduce errors in the interpolated value. */
2677 p->clut->re_set_rspl(
2678 p->clut, /* this */
2679 0, /* Combination of flags */
2680 (void *)p, /* Opaque function context */
2681 conv_rspl_out /* Function to set from */
2682 );
2683
2684 /* Fix absolute conversions to leave absolute response unchanged. */
2685 icmAry2XYZ(_wp, wp); /* Actual white point */
2686 if (p->picc != NULL) {
2687 p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, p->toAbs, p->fromAbs, icmD50, _wp);
2688 } else {
2689 icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
2690 icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
2691 }
2692
2693 if (p->verb) {
2694 double labwp[3];
2695
2696 /* Lookup white again */
2697 xfit_inpscurves(p, wcc.p, dw);
2698 p->clut->interp(p->clut, &wcc);
2699 xfit_outcurves(p, wcc.v, wcc.v);
2700 if (p->flags & XFIT_OUT_LAB)
2701 icmLab2XYZ(&icmD50, wcc.v, wcc.v);
2702 icmXYZ2Lab(&icmD50, labwp, wcc.v);
2703 printf("After fine tune, rel WP = XYZ %f %f %f, Lab %f %f %f\n",
2704 wcc.v[0], wcc.v[1], wcc.v[2], labwp[0], labwp[1], labwp[2]);
2705 printf(" abs WP = XYZ %s, Lab %s\n", icmPdv(3, wp), icmPLab(wp));
2706 }
2707 }
2708
2709 /* Create default wpscale */
2710 if (wpscale < 0.0) {
2711 wpscale = 1.0;
2712 } else {
2713 if (p->verb) {
2714 printf("White manual point scale %f\n", wpscale);
2715 }
2716 }
2717
2718 /* If we are going to auto scale the WP to avoid clipping */
2719 /* cLUT values above the WP: */
2720 if ((p->flags & XFIT_OUT_WP_REL_US) == XFIT_OUT_WP_REL_US) {
2721 co wcc;
2722 double bw[3];
2723 icmXYZNumber _wp;
2724 double uswpscale = 1.0;
2725 double mxd, mxY;
2726 double ndw[3];
2727
2728 /* See what device space gamut boundary white (ie. 1,1,1) maps to */
2729 xfit_inpscurves(p, wcc.p, dgw);
2730 p->clut->interp(p->clut, &wcc);
2731 xfit_outcurves(p, wcc.v, wcc.v);
2732 if (p->flags & XFIT_OUT_LAB)
2733 icmLab2XYZ(&icmD50, wcc.v, wcc.v);
2734 icmMulBy3x3(wcc.v, p->toAbs, wcc.v); /* Convert to absolute */
2735
2736 mxY = wcc.v[1];
2737 icmCpy3(bw, wcc.v);
2738 //printf("~1 1,1,1 Y = %f\n",wcc.v[1]);
2739
2740 /* See what the device white point value scaled to 1 produces */
2741 mxd = -1.0;
2742 for (e = 0; e < p->di; e++) {
2743 if (dw[e] > mxd)
2744 mxd = dw[e];
2745 }
2746 for (e = 0; e < p->di; e++)
2747 ndw[e] = dw[e]/mxd;
2748
2749 xfit_inpscurves(p, wcc.p, ndw);
2750 p->clut->interp(p->clut, &wcc);
2751 xfit_outcurves(p, wcc.v, wcc.v);
2752 if (p->flags & XFIT_OUT_LAB)
2753 icmLab2XYZ(&icmD50, wcc.v, wcc.v);
2754 icmMulBy3x3(wcc.v, p->toAbs, wcc.v); /* Convert to absolute */
2755
2756 //printf("~1 ndw = %f %f %f Y = %f\n",ndw[0],ndw[1],ndw[2],wcc.v[1]);
2757 if (wcc.v[1] > mxY) {
2758 mxY = wcc.v[1];
2759 icmCpy3(bw, wcc.v);
2760 }
2761
2762 /* Compute WP scale factor needed to fit mxY */
2763 if (mxY > wp[1]) {
2764 uswpscale = mxY/wp[1];
2765 wpscale *= uswpscale;
2766 if (p->verb) {
2767 printf("Dev boundary white XYZ %s, scale WP by %f, total WP scale %f\n",
2768 icmPdv(3, bw), uswpscale, wpscale);
2769 }
2770 }
2771 }
2772
2773 /* If the scaled WP would have Y > 1.0, clip it to 1.0 */
2774 if (p->flags & XFIT_CLIP_WP) {
2775
2776 if ((wp[1] * wpscale) > 1.0) {
2777 wpscale = 1.0/wp[1]; /* Make wp Y = 1.0 */
2778 if (p->verb) {
2779 printf("WP Y would ve > 1.0. scale by %f to clip it\n",wpscale);
2780 }
2781 }
2782 }
2783
2784 /* Apply our total wp scale factor */
2785 if (wpscale != 1.0) {
2786 icmXYZNumber _wp;
2787
2788 /* Create inverse scaling matrix for relative rspl data */
2789 icmSetUnity3x3(p->cmat);
2790 icmScale3x3(p->cmat, p->cmat, 1.0/wpscale);
2791
2792 /* Inverse scale the rspl */
2793 p->clut->re_set_rspl(
2794 p->clut, /* this */
2795 0, /* Combination of flags */
2796 (void *)p, /* Opaque function context */
2797 conv_rspl_out /* Function to set from */
2798 );
2799
2800 /* Scale the WP */
2801 icmScale3(wp, wp, wpscale);
2802
2803 /* return scaled white point to caller */
2804 icmAry2Ary(p->wp, wp);
2805
2806 /* Fix absolute conversions to leave absolute response unchanged. */
2807 icmAry2XYZ(_wp, wp); /* Actual white point */
2808 icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
2809 icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
2810 }
2811
2812 /* Clip any values in the grid over D50 L to D50 */
2813 if ((p->flags & XFIT_OUT_WP_REL_C) == XFIT_OUT_WP_REL_C) {
2814
2815 /* Compute the rspl D50 value to avoid calc in clip_rspl_out() */
2816 if (p->flags & XFIT_OUT_LAB) {
2817 p->cmat[0][0] = 100.0;
2818 p->cmat[0][1] = 0.0;
2819 p->cmat[0][2] = 0.0;
2820 } else {
2821 icmXYZ2Ary(p->cmat[0], icmD50);
2822 }
2823 xfit_invoutcurves(p, p->cmat[0], p->cmat[0]);
2824
2825 if (p->verb)
2826 printf("Clipping any cLUT grid points with Y > 1 to D50\n");
2827
2828 p->clut->re_set_rspl(
2829 p->clut, /* this */
2830 0, /* Combination of flags */
2831 (void *)p, /* Opaque function context */
2832 clip_rspl_out /* Function to set from */
2833 );
2834 }
2835
2836 /* Force black point to given value */
2837 // if (bpo != NULL) {
2838 // co tv;
2839 // int rv;
2840 //
2841 // xfit_inpscurves(p, tv.p, bpo->p);
2842 //
2843 // xfit_XYZ_abs_to_rel(p, tv.v, bpo->v);
2844 // xfit_invoutcurves(p, tv.v, tv.v);
2845 //printf("~1 xfit: fine after curves black at %f %f %f to %f %f %f\n",
2846 //tv.p[0], tv.p[1], tv.p[2], tv.v[0], tv.v[1], tv.v[2]);
2847 // rv = p->clut->tune_value(p->clut, &tv);
2848 // if (rv != 0)
2849 // warning("Black Point Override failed - clipping");
2850 // }
2851
2852 #ifdef SPECIAL_FORCE
2853 /* Replace the rspl nodes with ones directly computed */
2854 /* from the synthetic linear RGB->XYZ model */
2855 {
2856 double twp[3];
2857 icmXYZNumber _wp;
2858
2859 /* See what current device white maps to */
2860 twp[0] = twp[1] = twp[2] = 1.0;
2861 domodel(twp, twp);
2862
2863 icmAry2XYZ(_wp, twp);
2864
2865 /* Matrix needed to correct to D50 */
2866 icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->cmat);
2867
2868 p->clut->re_set_rspl(
2869 p->clut, /* this */
2870 0, /* Combination of flags */
2871 (void *)p, /* Opaque function context */
2872 set_rspl_out1 /* Function to set from */
2873 );
2874 }
2875 #endif /* SPECIAL_FORCE */
2876
2877 /* Evaluate the residual error now, with the rspl in place */
2878 #ifdef DEBUG_PLOT
2879 {
2880 int ee;
2881 double maxe = 0.0;
2882 double avee = 0.0;
2883
2884 /* Allocate in->rout duplicate point set */
2885 if (p->rpoints == NULL) {
2886 if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
2887 return 1;
2888 }
2889
2890 /* Create a set of 1D rspl setup points that contains */
2891 /* the residual error */
2892 for (i = 0; i < p->nodp; i++) {
2893 double tv[MXDO]; /* Target output value */
2894 double mv[MXDO]; /* Model output value */
2895 co pp;
2896 double ev;
2897
2898 xfit_abs_to_rel(p, tv, p->ipoints[i].v);
2899
2900 for (e = 0; e < p->di; e++)
2901 pp.p[e] = p->ipoints[i].p[e];
2902
2903 xfit_inpscurves(p, pp.p, pp.p);
2904 p->clut->interp(p->clut, &pp);
2905 xfit_outcurves(p, pp.v, pp.v);
2906
2907 for (f = 0; f < p->fdi; f++)
2908 mv[f] = pp.v[f];
2909
2910 /* Evaluate the residual error suqared */
2911 if (p->flags & XFIT_FM_INPUT) {
2912 double pp[MXDI];
2913 for (e = 0; e < di; e++)
2914 pp[e] = p->ipoints[i].p[e];
2915 for (f = 0; f < fdi; f++) {
2916 double t1 = tv[f] - mv[f]; /* Error in output */
2917
2918 /* Create input point offset by equivalent delta to output point */
2919 /* error, in proportion to the partial derivatives for that output. */
2920 for (e = 0; e < di; e++)
2921 pp[e] += t1 * p->piv[i].ide[f][e];
2922 }
2923 ev = p->to_de2(p->cntx2, pp, p->ipoints[i].p);
2924 } else {
2925 ev = p->to_de2(p->cntx2, mv, tv);
2926 }
2927 //printf("~1 point %d, loc %f %f %f %f\n",i, p->rpoints[i].p[0], p->rpoints[i].p[1], p->rpoints[i].p[2], p->rpoints[i].p[3]);
2928 //printf(" targ %f %f %f, is %f %f %f\n", tv[0], tv[1], tv[2], mv[0], mv[1], mv[2]);
2929
2930 p->rpoints[i].v[0] = ev;
2931 p->rpoints[i].w = p->ipoints[i].w;
2932
2933 ev = sqrt(ev);
2934 if (ev > maxe)
2935 maxe = ev;
2936 avee += ev;
2937 }
2938 printf("Max resid err = %f, avg err = %f\n",maxe, avee/(double)p->nodp);
2939
2940 /* Evaluate each input axis in turn */
2941 for (ee = 0; ee < p->di; ee++) {
2942 rspl *resid;
2943 double imin[1],imax[1],omin[1],omax[1];
2944 int resres[1] = { 1024 };
2945
2946 /* Create a rspl that gives the residual error squared */
2947 /* vs. the axis value */
2948 for (i = 0; i < p->nodp; i++)
2949 p->rpoints[i].p[0] = p->ipoints[i].p[ee];
2950
2951 imin[0] = in_min[ee];
2952 imax[0] = in_max[ee];
2953 omin[0] = 0.0;
2954 omax[0] = 0.0;
2955
2956 if ((resid = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL)
2957 return 1;
2958
2959 resid->fit_rspl_w(resid, RSPLFLAGS, p->rpoints, p->nodp, imin, imax, resres,
2960 omin, omax, 2.0, NULL, NULL);
2961
2962 {
2963 #define XRES 100
2964 double xx[XRES];
2965 double y1[XRES];
2966
2967 printf("Finale residual error vs. input channel %d\n",ee);
2968 for (i = 0; i < XRES; i++) {
2969 co pp;
2970 double x;
2971 x = i/(double)(XRES-1);
2972 xx[i] = x = x * (imax[0] - imin[0]) + imin[0];
2973 pp.p[0] = xx[i];
2974 resid->interp(resid, &pp);
2975 y1[i] = sqrt(fabs(pp.v[0]));
2976 }
2977 do_plot(xx,y1,NULL,NULL,XRES);
2978 }
2979 resid->del(resid);
2980 }
2981 }
2982 #endif /* DEBUG_PLOT */
2983
2984 /* Special test code to figure out what's wrong with position curves */
2985 #ifdef NEVER
2986 {
2987 double rgb[3], xyz[3];
2988 co pp;
2989
2990 extern int rspldb;
2991 db = 1;
2992 rspldb = 1;
2993 printf("~1 gres = %d %d %d\n", gres[0], gres[1], gres[2]);
2994 for (i = 0; i < 3; i++) {
2995
2996 if (i == 0) {
2997 /* Test point on a grid point */
2998 printf("\n~1 ##########################################\n");
2999 printf("~1 testing input at first diagonal grid point\n");
3000
3001 rgb[0] = 1.0/(gres[0]-1.0);
3002 rgb[1] = 1.0/(gres[0]-1.0);
3003 rgb[2] = 1.0/(gres[0]-1.0);
3004
3005 printf("~1 target rgb' = %f %f %f\n", rgb[0], rgb[1], rgb[2]);
3006 xfit_invinpscurves(p, rgb, rgb);
3007
3008 } else if (i == 1) {
3009 /* Test point half way through a grid point */
3010 printf("\n~1 ##########################################\n");
3011 printf("\n~1 testing half way through diagonal grid point\n");
3012
3013 rgb[0] = 0.5/(gres[0]-1.0);
3014 rgb[1] = 0.5/(gres[0]-1.0);
3015 rgb[2] = 0.5/(gres[0]-1.0);
3016
3017 printf("~1 target rgb' = %f %f %f\n", rgb[0], rgb[1], rgb[2]);
3018 xfit_invinpscurves(p, rgb, rgb);
3019
3020 } else {
3021 printf("\n~1 ##########################################\n");
3022 printf("\n~1 testing worst case point\n");
3023
3024 rgb[0] = 0.039915;
3025 rgb[1] = 0.053148;
3026 rgb[2] = 0.230610;
3027 }
3028
3029 pp.p[0] = rgb[0];
3030 pp.p[1] = rgb[1];
3031 pp.p[2] = rgb[2];
3032
3033 printf("~1 rgb = %f %f %f\n", pp.p[0], pp.p[1], pp.p[2]);
3034
3035 xfit_inpscurves(p, pp.p, pp.p);
3036
3037 printf("~1 rgb' = %f %f %f\n", pp.p[0], pp.p[1], pp.p[2]);
3038
3039 p->clut->interp(p->clut, &pp);
3040
3041 printf("~1 xyz' = %f %f %f\n", pp.v[0], pp.v[1], pp.v[2]);
3042
3043 xfit_outcurves(p, pp.v, pp.v);
3044
3045 printf("~1 xyz = %f %f %f\n", pp.v[0], pp.v[1], pp.v[2]);
3046
3047 /* Synthetic linear rgb->XYZ model */
3048 domodel(xyz, rgb);
3049
3050 /* Apply abs->rel white point adjustment */
3051 icmMulBy3x3(xyz, p->mat, xyz);
3052
3053 printf("~1 ref = %f %f %f, de = %f\n", xyz[0], xyz[1], xyz[2],icmXYZLabDE(&icmD50,xyz,pp.v));
3054 }
3055
3056 db = 0;
3057 rspldb = 0;
3058 // exit(0);
3059 }
3060 #endif /* NEVER */
3061
3062 }
3063 return 0;
3064 }
3065
3066 /* We're done with an xfit */
xfit_del(xfit * p)3067 static void xfit_del(xfit *p) {
3068 if (p->v != NULL)
3069 free(p->v);
3070 if (p->wv != NULL)
3071 free(p->wv);
3072 if (p->sa != NULL)
3073 free(p->sa);
3074 if (p->rpoints != NULL)
3075 free(p->rpoints);
3076 if (p->piv != NULL)
3077 free(p->piv);
3078 if (p->uerrv != NULL)
3079 free(p->uerrv);
3080 free(p);
3081 }
3082
3083 /* Create a transform fitting object */
3084 /* return NULL on error */
new_xfit(icc * picc)3085 xfit *new_xfit(
3086 icc *picc /* ICC profile used to set cone space matrix, NULL for Bradford. */
3087 ) {
3088 xfit *p;
3089
3090 if ((p = (xfit *)calloc(1, sizeof(xfit))) == NULL) {
3091 return NULL;
3092 }
3093
3094 p->picc = picc;
3095
3096 /* Set method pointers */
3097 p->fit = xfit_fit;
3098 p->incurve = xfit_inpscurve;
3099 p->invincurve = xfit_invinpscurve;
3100 p->outcurve = xfit_outcurve;
3101 p->invoutcurve = xfit_invoutcurve;
3102 p->del = xfit_del;
3103
3104 return p;
3105 }
3106
3107
3108
3109