1
2 /*
3 * Model Printer Profile Lookup test utility.
4 *
5 * Author: Graeme W. Gill
6 * Date: 2002/12/30
7 * Version: 1.00
8 *
9 * Copyright 2002, 2003 Graeme W. Gill
10 * All rights reserved.
11 * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
12 * see the License.txt file for licencing details.
13 *
14 */
15
16 /* TTBD:
17 *
18 */
19
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <stdarg.h>
23 #include <fcntl.h>
24 #include <string.h>
25 #include <math.h>
26 #include "numlib.h"
27 #include "xicc.h"
28 #include "counters.h"
29 #include "vrml.h"
30 #include "ui.h"
31
usage(void)32 void usage(void) {
33 fprintf(stderr,"Translate colors through an MPP profile, V1.00\n");
34 fprintf(stderr,"Author: Graeme W. Gill\n");
35 fprintf(stderr,"usage: mpplu [-v] [-f func] [-i intent] [-o order] profile\n");
36 fprintf(stderr," -v Verbose\n");
37 fprintf(stderr," -f function f = forward, b = backwards\n");
38 fprintf(stderr," -p oride x = XYZ_PCS, l = Lab_PCS, y = Yxy, s = spectral,\n");
39 fprintf(stderr," -l limit override default ink limit, 1 - N00%%\n");
40 fprintf(stderr," -i illum Choose illuminant for print/transparency spectral data:\n");
41 fprintf(stderr," A, C, D50 (def.), D50M2, D65, F5, F8, F10 or file.sp\n");
42 fprintf(stderr," -o observ Choose CIE Observer for spectral data:\n");
43 fprintf(stderr," 1931_2 (def), 1964_10, S&B 1955_2, shaw, J&V 1978_2\n");
44 fprintf(stderr," -u Use Fluorescent Whitening Agent compensation\n");
45 fprintf(stderr," -g Create gamut output\n");
46 fprintf(stderr," -w Create gamut %s as well\n",vrml_format());
47 fprintf(stderr," -n Don't add %s axes\n",vrml_format());
48 fprintf(stderr," -a n Gamut transparency level\n");
49 fprintf(stderr," -d n Gamut surface detail level\n");
50 fprintf(stderr," -t num Invoke debugging test code \"num\" 1..n\n");
51 fprintf(stderr," 1 - check partial derivative for device input\n");
52 fprintf(stderr," 2 - create overlap diagnostic %s gamut surface\n",vrml_format());
53 fprintf(stderr,"\n");
54 fprintf(stderr," The colors to be translated should be fed into stdin,\n");
55 fprintf(stderr," one input color per line, white space separated.\n");
56 fprintf(stderr," A line starting with a # will be ignored.\n");
57 fprintf(stderr," A line not starting with a number will terminate the program.\n");
58 exit(1);
59 }
60
61 static void diag_gamut(mpp *p, double gamres, int doaxes, double trans, char *outname);
62 static void mpp_rev(mpp *mppo, double limit, double *out, double *in);
63
64 int
main(int argc,char * argv[])65 main(int argc, char *argv[]) {
66 int fa,nfa; /* argument we're looking at */
67 char prof_name[100];
68 mpp *mppo;
69 int verb = 0;
70 int test = 0; /* special test code */
71 int dogam = 0; /* Create gamut */
72 int dowrl = 0; /* Create VRML/X3D gamut */
73 int doaxes = 1; /* Create VRML/X3D axes */
74 double trans = 0.0; /* Transparency */
75 double gamres = 0.0; /* Gamut resolution */
76 int repYxy = 0; /* Report Yxy */
77 int repSpec = 0; /* Report Spectral */
78 int bwd = 0; /* Do reverse lookup */
79 double dlimit; /* Device ink limit */
80 double limit = -1.0; /* Used ink limit */
81 int display = 0; /* NZ if display type */
82 int spec = 0; /* Use spectral data flag */
83 int spec_n; /* Number of spectral bands, 0 if not valid */
84 double spec_wl_short; /* First reading wavelength in nm (shortest) */
85 double spec_wl_long; /* Last reading wavelength in nm (longest) */
86 int fwacomp = 0; /* FWA compensation */
87 icxIllumeType illum = icxIT_default; /* Spectral defaults */
88 xspect cust_illum; /* Custom illumination spectrum */
89 icxObserverType observ = icxOT_default;
90 char buf[200];
91 double in[MAX_CHAN], out[MAX_CHAN];
92 int rv = 0;
93
94 inkmask imask; /* Device Ink mask */
95 char *ident = NULL; /* Device colorspec description */
96 icColorSpaceSignature pcss; /* Type of PCS space */
97 int devn, pcsn; /* Number of components */
98
99 icColorSpaceSignature pcsor = icSigLabData; /* Default */
100
101 error_program = argv[0];
102
103 if (argc < 2)
104 usage();
105
106 /* Process the arguments */
107 for(fa = 1;fa < argc;fa++) {
108 nfa = fa; /* skip to nfa if next argument is used */
109 if (argv[fa][0] == '-') { /* Look for any flags */
110 char *na = NULL; /* next argument after flag, null if none */
111
112 if (argv[fa][2] != '\000')
113 na = &argv[fa][2]; /* next is directly after flag */
114 else {
115 if ((fa+1) < argc) {
116 if (argv[fa+1][0] != '-') {
117 nfa = fa + 1;
118 na = argv[nfa]; /* next is seperate non-flag argument */
119 }
120 }
121 }
122
123 if (argv[fa][1] == '?')
124 usage();
125
126 /* Verbosity */
127 else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') {
128 verb = 1;
129 }
130 /* function */
131 else if (argv[fa][1] == 'f' || argv[fa][1] == 'F') {
132 fa = nfa;
133 if (na == NULL) usage();
134 switch (na[0]) {
135 case 'f':
136 case 'F':
137 bwd = 0;
138 break;
139 case 'b':
140 case 'B':
141 bwd = 1;
142 break;
143 default:
144 usage();
145 }
146 }
147
148 /* PCS override */
149 else if (argv[fa][1] == 'p' || argv[fa][1] == 'P') {
150 fa = nfa;
151 if (na == NULL) usage();
152 switch (na[0]) {
153 case 'x':
154 case 'X':
155 pcsor = icSigXYZData;
156 repYxy = repSpec = 0;
157 break;
158 case 'l':
159 case 'L':
160 pcsor = icSigLabData;
161 repYxy = repSpec = 0;
162 break;
163 case 'y':
164 case 'Y':
165 pcsor = icSigXYZData;
166 repYxy = 1;
167 repSpec = 0;
168 break;
169 case 's':
170 case 'S':
171 pcsor = icSigXYZData;
172 repYxy = 0;
173 repSpec = 1;
174 spec = 1;
175 break;
176 default:
177 usage();
178 }
179 }
180
181 /* Ink Limit */
182 else if (argv[fa][1] == 'l' || argv[fa][1] == 'L') {
183 fa = nfa;
184 if (na == NULL) usage();
185 limit = atof(na);
186 }
187
188 /* Spectral Illuminant type */
189 else if (argv[fa][1] == 'i' || argv[fa][1] == 'I') {
190 fa = nfa;
191 if (na == NULL) usage();
192 if (strcmp(na, "A") == 0) {
193 spec = 1;
194 illum = icxIT_A;
195 } else if (strcmp(na, "C") == 0) {
196 spec = 1;
197 illum = icxIT_C;
198 } else if (strcmp(na, "D50") == 0) {
199 spec = 1;
200 illum = icxIT_D50;
201 } else if (strcmp(na, "D50M2") == 0) {
202 spec = 1;
203 illum = icxIT_D50M2;
204 } else if (strcmp(na, "D65") == 0) {
205 spec = 1;
206 illum = icxIT_D65;
207 } else if (strcmp(na, "F5") == 0) {
208 spec = 1;
209 illum = icxIT_F5;
210 } else if (strcmp(na, "F8") == 0) {
211 spec = 1;
212 illum = icxIT_F8;
213 } else if (strcmp(na, "F10") == 0) {
214 spec = 1;
215 illum = icxIT_F10;
216 } else { /* Assume it's a filename */
217 inst_meas_type mt;
218
219 spec = 1;
220 illum = icxIT_custom;
221 if (read_xspect(&cust_illum, &mt, na) != 0)
222 usage();
223
224 if (mt != inst_mrt_none
225 && mt != inst_mrt_emission
226 && mt != inst_mrt_ambient
227 && mt != inst_mrt_emission_flash
228 && mt != inst_mrt_ambient_flash) {
229 error("Custom illuminant '%s' is wrong measurement type",na);
230 }
231 }
232 }
233
234 /* Spectral Observer type */
235 else if (argv[fa][1] == 'o' || argv[fa][1] == 'O') {
236 fa = nfa;
237 if (na == NULL) usage();
238 if (strcmp(na, "1931_2") == 0) { /* Classic 2 degree */
239 spec = 1;
240 observ = icxOT_CIE_1931_2;
241 } else if (strcmp(na, "1964_10") == 0) { /* Classic 10 degree */
242 spec = 1;
243 observ = icxOT_CIE_1964_10;
244 } else if (strcmp(na, "1955_2") == 0) { /* Stiles and Burch 1955 2 degree */
245 spec = 1;
246 observ = icxOT_Stiles_Burch_2;
247 } else if (strcmp(na, "1978_2") == 0) { /* Judd and Voss 1978 2 degree */
248 spec = 1;
249 observ = icxOT_Judd_Voss_2;
250 } else if (strcmp(na, "shaw") == 0) { /* Shaw and Fairchilds 1997 2 degree */
251 spec = 1;
252 observ = icxOT_Shaw_Fairchild_2;
253 } else
254 usage();
255 }
256
257 /* Fluerescent Whitner compensation */
258 else if (argv[fa][1] == 'u' || argv[fa][1] == 'U')
259 fwacomp = 1;
260
261 /* Gamut plot */
262 else if (argv[fa][1] == 'g' || argv[fa][1] == 'G')
263 dogam = 1;
264
265 /* VRML/X3D plot */
266 else if (argv[fa][1] == 'w' || argv[fa][1] == 'W') {
267 dogam = 1;
268 dowrl = 1;
269 }
270
271 /* No VRML/X3D axes */
272 else if (argv[fa][1] == 'n' || argv[fa][1] == 'N') {
273 doaxes = 0;
274 }
275
276 /* Transparency */
277 else if (argv[fa][1] == 'a' || argv[fa][1] == 'A') {
278 fa = nfa;
279 if (na == NULL) usage();
280 trans = atof(na);
281 }
282
283 /* Surface Detail */
284 else if (argv[fa][1] == 'd' || argv[fa][1] == 'D') {
285 fa = nfa;
286 if (na == NULL) usage();
287 gamres = atof(na);
288 dogam = 1;
289 }
290
291 /* Test code */
292 else if (argv[fa][1] == 't' || argv[fa][1] == 'T') {
293 fa = nfa;
294 if (na == NULL) usage();
295 test = atoi(na);
296 }
297
298 else
299 usage();
300 } else
301 break;
302 }
303
304 if (fa >= argc || argv[fa][0] == '-') usage();
305 strcpy(prof_name,argv[fa]);
306
307 if ((mppo = new_mpp()) == NULL)
308 error ("Creation of MPP object failed");
309
310 if ((rv = mppo->read_mpp(mppo,prof_name)) != 0)
311 error ("%d, %s",rv,mppo->err);
312
313 mppo->get_info(mppo, &imask, &devn, &dlimit, &spec_n, &spec_wl_short, &spec_wl_long, NULL, &display);
314 ident = icx_inkmask2char(imask, 1);
315
316 if (verb) {
317 printf("MPP profile with %d colorants, type %s, TAC %f\n",devn,ident, dlimit);
318 if (display)
319 printf("MPP profile is for a display type device\n");
320 }
321
322 if (limit <= 0.0 || dlimit < limit)
323 limit = dlimit;
324
325 pcss = pcsor;
326 pcsn = 3;
327
328 if (spec && spec_n == 0) {
329 error("Spectral profile needed for spectral result, custom illuminant, observer or FWA");
330 }
331
332 /* Select CIE return value details */
333 if ((rv = mppo->set_ilob(mppo, illum, &cust_illum, observ, NULL, pcss, fwacomp)) != 0) {
334 if (rv == 1)
335 error("Spectral profile needed for custom illuminant, observer or FWA");
336 error("Error setting illuminant, observer, or FWA");
337 }
338
339 if (test != 0) {
340 printf("!!!!! Running special test code no %d !!!!!\n",test);
341
342 if (test == 1) {
343 double **dv, **rdv;
344
345 dv = dmatrix(0, pcsn-1, 0, devn-1);
346 rdv = dmatrix(0, pcsn-1, 0, devn-1);
347
348 printf("Checking partial derivative at each input value\n");
349 /* Process colors to translate */
350 for (;;) {
351 int i,j;
352 char *bp, *nbp;
353 double tout[MAX_CHAN];
354
355 /* Read in the next line */
356 if (fgets(buf, 200, stdin) == NULL)
357 break;
358 if (buf[0] == '#') {
359 fprintf(stdout,"%s\n",buf);
360 continue;
361 }
362 /* For each input number */
363 for (nbp = buf, i = 0; i < MAX_CHAN; i++) {
364 bp = nbp;
365 in[i] = strtod(bp, &nbp);
366 if (nbp == bp)
367 break; /* Failed */
368 }
369 if (i == 0)
370 break;
371
372 /* Do conversion */
373 mppo->lookup(mppo, out, in);
374 mppo->dlookup(mppo, out, dv, in);
375
376 for (j = 0; j < devn; j++) {
377 double del = 1e-9;
378
379 if (in[j] > 0.5)
380 del = -del;
381
382 in[j] += del;
383 mppo->lookup(mppo, tout, in);
384 in[j] -= del;
385
386 for (i = 0; i < pcsn; i++) {
387 rdv[i][j] = (tout[i] - out[i])/del;
388 }
389 }
390
391 /* Output the results */
392 for (j = 0; j < devn; j++) {
393 if (j > 0)
394 fprintf(stdout," %f",in[j]);
395 else
396 fprintf(stdout,"%f",in[j]);
397 }
398 printf(" [%s] -> ", ident);
399
400 for (j = 0; j < pcsn; j++) {
401 if (j > 0)
402 fprintf(stdout," %f",out[j]);
403 else
404 fprintf(stdout,"%f",out[j]);
405 }
406 printf(" [%s]\n", icm2str(icmColorSpaceSignature, pcss));
407
408 /* Print the derivatives */
409 for (i = 0; i < pcsn; i++) {
410
411 printf("Output chan %d: ",i);
412 for (j = 0; j < devn; j++) {
413 if (j < (devn-1))
414 fprintf(stdout,"%f ref %f, ",dv[i][j], rdv[i][j]);
415 else
416 fprintf(stdout,"%f ref %f\n",dv[i][j], rdv[i][j]);
417 }
418 }
419 }
420
421 free_dmatrix(dv, 0, pcsn-1, 0, devn-1);
422 free_dmatrix(rdv, 0, pcsn-1, 0, devn-1);
423
424 } else if (test == 2) {
425 char *xl, gam_name[100];
426
427 strcpy(gam_name, prof_name);
428 if ((xl = strrchr(gam_name, '.')) == NULL) /* Figure where extention is */
429 xl = gam_name + strlen(gam_name);
430 xl[0] = '\000'; /* Remove extension */
431 diag_gamut(mppo, gamres, doaxes, trans, gam_name);
432
433 } else {
434 printf("Unknown test!\n");
435 }
436
437 } else if (dogam) {
438 gamut *gam;
439 char *xl, gam_name[100];
440 int docusps = 1;
441
442 if ((gam = mppo->get_gamut(mppo, gamres)) == NULL)
443 error("get_gamut failed\n");
444
445 strcpy(gam_name, prof_name);
446 if ((xl = strrchr(gam_name, '.')) == NULL) /* Figure where extention is */
447 xl = gam_name + strlen(gam_name);
448 strcpy(xl,".gam");
449 if (gam->write_gam(gam, gam_name))
450 error ("write gamut failed on '%s'",gam_name);
451
452 if (dowrl) {
453 xl[0] = '\000';
454 if (gam->write_vrml(gam,gam_name, doaxes, docusps))
455 error ("write vrml failed on '%s%s'",gam_name,vrml_ext());
456 }
457
458 gam->del(gam);
459
460 } else {
461 /* Normal color lookup */
462
463 if (repYxy) { /* report Yxy rather than XYZ */
464 if (pcss == icSigXYZData)
465 pcss = icSigYxyData;
466 }
467
468 /* Process colors to translate */
469 for (;;) {
470 int i,j;
471 char *bp, *nbp;
472
473 /* Read in the next line */
474 if (fgets(buf, 200, stdin) == NULL)
475 break;
476 if (buf[0] == '#') {
477 fprintf(stdout,"%s\n",buf);
478 continue;
479 }
480 /* For each input number */
481 for (nbp = buf, i = 0; i < MAX_CHAN; i++) {
482 bp = nbp;
483 in[i] = strtod(bp, &nbp);
484 if (nbp == bp)
485 break; /* Failed */
486 }
487 if (i == 0)
488 break;
489
490 if (!bwd) {
491
492 if (repSpec) {
493 xspect ospec;
494
495 /* Do lookup of spectrum */
496 mppo->lookup_spec(mppo, &ospec, in);
497
498 /* Output the results */
499 for (j = 0; j < devn; j++) {
500 if (j > 0)
501 fprintf(stdout," %f",in[j]);
502 else
503 fprintf(stdout,"%f",in[j]);
504 }
505 printf(" [%s] -> ", ident);
506
507 for (j = 0; j < spec_n; j++) {
508 if (j > 0)
509 fprintf(stdout," %f",ospec.spec[j]);
510 else
511 fprintf(stdout,"%f",ospec.spec[j]);
512 }
513
514 printf(" [%3.0f .. %3.0f nm]\n", spec_wl_short, spec_wl_long);
515
516 } else {
517
518 /* Do conversion */
519 mppo->lookup(mppo, out, in);
520
521 /* Output the results */
522 for (j = 0; j < devn; j++) {
523 if (j > 0)
524 fprintf(stdout," %f",in[j]);
525 else
526 fprintf(stdout,"%f",in[j]);
527 }
528 printf(" [%s] -> ", ident);
529
530 if (repYxy && pcss == icSigYxyData) {
531 double X = out[0];
532 double Y = out[1];
533 double Z = out[2];
534 double sum = X + Y + Z;
535 if (sum < 1e-6) {
536 out[0] = out[1] = out[2] = 0.0;
537 } else {
538 out[0] = Y;
539 out[1] = X/sum;
540 out[2] = Y/sum;
541 }
542 }
543 for (j = 0; j < pcsn; j++) {
544 if (j > 0)
545 fprintf(stdout," %f",out[j]);
546 else
547 fprintf(stdout,"%f",out[j]);
548 }
549
550 printf(" [%s]\n", icm2str(icmColorSpaceSignature, pcss));
551 }
552
553 } else { /* Do a reverse lookup */
554
555 if (repYxy && pcss == icSigYxyData) {
556 double Y = in[0];
557 double x = in[1];
558 double y = in[2];
559 double z = 1.0 - x - y;
560 double sum;
561 if (y < 1e-6) {
562 in[0] = in[1] = in[2] = 0.0;
563 } else {
564 sum = Y/y;
565 in[0] = x * sum;
566 in[1] = Y;
567 in[2] = z * sum;
568 }
569 }
570
571 /* Do conversion */
572 mpp_rev(mppo, limit, out, in);
573
574 /* Output the results */
575 for (j = 0; j < pcsn; j++) {
576 if (j > 0)
577 fprintf(stdout," %f",in[j]);
578 else
579 fprintf(stdout,"%f",in[j]);
580 }
581 printf(" [%s] -> ", icm2str(icmColorSpaceSignature, pcss));
582
583 for (j = 0; j < devn; j++) {
584 if (j > 0)
585 fprintf(stdout," %f",out[j]);
586 else
587 fprintf(stdout,"%f",out[j]);
588 }
589
590 printf(" [%s]\n", ident);
591 }
592 }
593 }
594
595 free(ident);
596 mppo->del(mppo);
597
598 return 0;
599 }
600
601
602 /* -------------------------------------------- */
603 /* Code for special gamut surface plot */
604
605 #define GAMUT_LCENT 50
606
607 /* Create a diagnostic gamut, illustrating */
608 /* device space "fold-over" */
609 /* (This will be in the current PCS, but assumed to be Lab) */
diag_gamut(mpp * p,double detail,int doaxes,double trans,char * outname)610 static void diag_gamut(
611 mpp *p, /* This */
612 double detail, /* Gamut resolution detail */
613 int doaxes,
614 double trans, /* Transparency */
615 char *outname /* Output VRML/X3D file (no extension) */
616 ) {
617 int i, j;
618 vrml *wrl;
619 int vix; /* Vertex index */
620 DCOUNT(coa, MAX_CHAN, p->n, 0, 0, 2);
621 double col[MPP_MXCCOMB][3]; /* Color asigned to each major vertex */
622 int res;
623
624 /* Asign some colors to the combination nodes */
625 for (i = 0; i < p->nn; i++) {
626 int a, b, c, j;
627 double h;
628
629 j = (i ^ 0x5a5a5a5a) % p->nn;
630 h = (double)j/(p->nn-1);
631
632 /* Make fully saturated with chosen hue */
633 if (h < 1.0/3.0) {
634 a = 0;
635 b = 1;
636 c = 2;
637 } else if (h < 2.0/3.0) {
638 a = 1;
639 b = 2;
640 c = 0;
641 h -= 1.0/3.0;
642 } else {
643 a = 2;
644 b = 0;
645 c = 1;
646 h -= 2.0/3.0;
647 }
648 h *= 3.0;
649
650 col[i][a] = (1.0 - h);
651 col[i][b] = h;
652 col[i][c] = d_rand(0.0, 1.0);
653 }
654
655 if (detail > 0.0)
656 res = (int)(100.0/detail); /* Establish an appropriate sampling density */
657 else
658 res = 4;
659
660 if (res < 2)
661 res = 2;
662
663 if ((wrl = new_vrml(outname, doaxes, vrml_lab)) == NULL)
664 error("new_vrml faile for file '%s%s'",outname,vrml_ext());
665
666 wrl->start_line_set(wrl, 0);
667
668 /* Itterate over all the faces in the device space */
669 /* generating the vertx positions. */
670 DC_INIT(coa);
671 vix = 0;
672 while(!DC_DONE(coa)) {
673 int e, m1, m2;
674 double in[MAX_CHAN];
675 double out[3];
676 double sum;
677 double xb, yb;
678
679 /* Scan only device surface */
680 for (m1 = 0; m1 < p->n; m1++) {
681 if (coa[m1] != 0)
682 continue;
683
684 for (m2 = m1 + 1; m2 < p->n; m2++) {
685 int x, y;
686
687 if (coa[m2] != 0)
688 continue;
689
690 for (sum = 0.0, e = 0; e < p->n; e++)
691 in[e] = (double)coa[e]; /* Base value */
692
693 /* Scan over 2D device space face */
694 for (x = 0; x < res; x++) { /* step over surface */
695 xb = in[m1] = x/(res - 1.0);
696 for (y = 0; y < res; y++) {
697 int v0, v1, v2, v3;
698 double rgb[3];
699 yb = in[m2] = y/(res - 1.0);
700
701 /* Lookup PCS value */
702 p->lookup(p, out, in);
703
704 /* Create a color */
705 for (v0 = 0, e = 0; e < p->n; e++)
706 v0 |= coa[e] ? (1 << e) : 0; /* Binary index */
707
708 v1 = v0 | (1 << m2); /* Y offset */
709 v2 = v0 | (1 << m2) | (1 << m1); /* X+Y offset */
710 v3 = v0 | (1 << m1); /* Y offset */
711
712 /* Linear interp between the main verticies */
713 for (j = 0; j < 3; j++) {
714 rgb[j] = (1.0 - yb) * (1.0 - xb) * col[v0][j]
715 + yb * (1.0 - xb) * col[v1][j]
716 + (1.0 - yb) * xb * col[v3][j]
717 + yb * xb * col[v2][j];
718 }
719
720 wrl->add_col_vertex(wrl, 0, out, rgb);
721 vix++;
722 }
723 }
724 }
725 }
726 /* Increment index within block */
727 DC_INC(coa);
728 }
729
730 /* Itterate over all the faces in the device space */
731 /* generating the quadrilateral indexes. */
732 DC_INIT(coa);
733 vix = 0;
734 while(!DC_DONE(coa)) {
735 int e, m1, m2;
736 double in[MAX_CHAN];
737 double sum;
738
739 /* Scan only device surface */
740 for (m1 = 0; m1 < p->n; m1++) {
741 if (coa[m1] != 0)
742 continue;
743
744 for (m2 = m1 + 1; m2 < p->n; m2++) {
745 int x, y;
746
747 if (coa[m2] != 0)
748 continue;
749
750 for (sum = 0.0, e = 0; e < p->n; e++)
751 in[e] = (double)coa[e]; /* Base value */
752
753 /* Scan over 2D device space face */
754 for (x = 0; x < res; x++) { /* step over surface */
755 for (y = 0; y < res; y++) {
756
757 if (x < (res-1) && y < (res-1)) {
758 int ix[4];
759 ix[0] = vix;
760 ix[1] = vix + 1;
761 ix[2] = vix + 1 + res;
762 ix[3] = vix + res;
763 wrl->add_quad(wrl, 0, ix);
764 }
765 vix++;
766 }
767 }
768 }
769 }
770 /* Increment index within block */
771 DC_INC(coa);
772 }
773
774 wrl->make_quads_vc(wrl, 0, trans);
775
776 if (wrl->flush(wrl) != 0)
777 error("Error closing output file '%s%s'",outname,vrml_ext());
778
779 wrl->del(wrl);
780 }
781
782 /* -------------------------------------------- */
783 /* Reverse lookup support */
784 /* This is for developing the appropriate reverse lookup */
785 /* code for xsep.c */
786
787 /*
788 * TTBD:
789 * Not handlink rule or separate black
790 * Not handling linearisation callback for ink limit.
791 * Not allowing for other possible secondary limits/goals/tunings.
792 */
793
794 /*
795 * Descriptionr:
796 *
797 * The primary reverse lookup optimisation goals are to remain within
798 * the device gamut, and to match the target PCS.
799 *
800 * The secondary optimisation goals for solving this under constrained
801 * problem can be chosen to a achieve a wide variety of possible aims.
802 *
803 * 1) One that is applicable to screened devices might be to use the extra
804 * inks to minimise the visiblility of screening. For screening resolutions
805 * above 50 lpi (equivalent to 100 dpi stocastic screens) only luminance
806 * contrast will be relavant, so priority to the inks closest to paper white
807 * measured purely by L distance is appropriate. (In practice I've used
808 * a measure that adds a small color distance component.)
809 *
810 * 2) Another possible goal would be to optimise fit to a desired spectral
811 * profile, if spectral information is avaiable as an aim. The goal would
812 * be best fit weighted to the visual sensitivity curve.
813 *
814 * 3) Another possible secondary goal might be to minimise the total
815 * amount of ink used, to minimise cost, drying time, media ink loading.
816 *
817 * 4) A fourth possible secondary goal might be to choose ink combinations
818 * that are the most robust given an error in device values.
819 *
820 * In practice these secondary goals need not be entirely exclusive,
821 * but the overall secondary goal could be a weighted blending between
822 * these goals. Overall the ink limit (TAC) is extremely important,
823 * since this will be the primary thing that stops large amounds of
824 * light ink being used at all times.
825 *
826 * Numerous other tweaks, limits or goals (such as secondary combination
827 * ink limits, exclusions such as Cyan/Oraange, Magenta/Green)
828 * could also be applied in the reverse optimisation routine.
829 *
830 */
831
832
833 #ifdef NEVER
834 /* These weights are the "theoreticaly correct" weightings, since */
835 /* at 50 lpi or higher, the color contrast sensitivity should be close to 0 */
836 # define L_WEIGHT 1.0
837 # define a_WEIGHT 0.0
838 # define b_WEIGHT 0.0
839 #else
840 /* These weights give us our "expected" ink ordering of */
841 /* Yellow, light Cyan/Magenta, Orange/Green, Cyan/Magenta, Black. */
842 # define L_WEIGHT 1.0
843 # define a_WEIGHT 0.4
844 # define b_WEIGHT 0.2
845 #endif
846
847 /* Start array entry */
848 typedef struct {
849 double dev[MAX_CHAN]; /* Device value */
850 double Lab[3]; /* Lab value */
851 double oerr; /* Order error */
852 } saent;
853
854 /* Context for reverse lookup */
855 typedef struct {
856 int pass;
857 int di; /* Number of device channels */
858 double Lab[3]; /* Lab target value */
859 void (*dev2lab) (mpp *d2lcntx, double *out, double *in); /* Device to Lab callback */
860 mpp *d2lcntx; /* dev2lab callback context */
861 double ilimit; /* Total limit */
862 int sord[MAX_CHAN]; /* Sorted order index */
863 double oweight[MAX_CHAN]; /* Order weighting (not used ?) */
864 } revlus;
865
866 /* Return the largest distance of the point outside the device gamut. */
867 /* This will be 0 if inside the gamut, and > 0 if outside. */
868 static double
dist2gamut(revlus * s,double * d)869 dist2gamut(revlus *s, double *d) {
870 int e;
871 int di = s->di;
872 double tt, dd = 0.0;
873 double ss = 0.0;
874
875 for (e = 0; e < di; e++) {
876 ss += d[e];
877
878 tt = 0.0 - d[e];
879 if (tt > 0.0) {
880 if (tt > dd)
881 dd = tt;
882 }
883 tt = d[e] - 1.0;
884 if (tt > 0.0) {
885 if (tt > dd)
886 dd = tt;
887 }
888 }
889 tt = ss - s->ilimit;
890 if (tt > 0.0) {
891 if (tt > dd)
892 dd = tt;
893 }
894 return dd;
895 }
896
897 /* Snap a point to the device gamut boundary. */
898 /* Return nz if it has been snapped. */
snap2gamut(revlus * s,double * d)899 static int snap2gamut(revlus *s, double *d) {
900 int e;
901 int di = s->di;
902 double dd; /* Smallest distance */
903 double ss; /* Sum */
904 int rv = 0;
905
906 /* Snap to ink limit first */
907 for (ss = 0.0, e = 0; e < di; e++)
908 ss += d[e];
909 dd = fabs(ss - s->ilimit);
910
911 if (dd < 0.0) {
912 int j;
913 for (j = 0; j < di; j++)
914 d[j] *= s->ilimit/ss; /* Snap to ink limit */
915 rv = 1;
916 }
917
918 /* Now snap to any other dimension */
919 for (e = 0; e < di; e++) {
920
921 dd = fabs(d[e] - 0.0);
922 if (dd < 0.0) {
923 d[e] = 0.0; /* Snap to orthogonal boundary */
924 rv = 1;
925 }
926 dd = fabs(1.0 - d[e]);
927 if (dd < 0.0) {
928 d[e] = 1.0; /* Snap to orthogonal boundary */
929 rv = 1;
930 }
931 }
932
933 return rv;
934 }
935
936 /* Reverse optimisation function handed to powell() */
revoptfunc(void * edata,double * v)937 static double revoptfunc(void *edata, double *v) {
938 revlus *s = (revlus *)edata;
939 double rv;
940
941 printf("~1 target %f %f %f\n",s->Lab[0],s->Lab[1],s->Lab[2]);
942
943 if ((rv = (dist2gamut(s, v))) > 0.0) {
944 // rv = rv * 1000.0 + 45000.0; /* Discourage being out of gamut */
945 rv = rv * 5e6; /* Discourage being out of gamut */
946
947 }
948 printf("~1 out of gamut error = %f\n",rv);
949 {
950 int j;
951 double oerr;
952 double Lab[3];
953 double tot;
954
955 /* Convert device to Lab */
956 s->dev2lab(s->d2lcntx, Lab, v);
957
958 /* Accumulate total delta E squared */
959 for (j = 0; j < 3; j++) {
960 double tt = s->Lab[j] - Lab[j];
961 rv += tt * tt;
962 }
963
964 printf("~1 Delta E squared = %f\n",rv);
965
966 /* Skip first 3 colorants */
967 oerr = tot = 0.0;
968 printf("oerr = %f\n",oerr);
969 for (j = 3; j < s->di; j++) {
970 int ix = s->sord[j]; /* Sorted order index */
971 double vv = v[ix];
972 double we = (double)j - 2.0; /* Increasing weight */
973
974 printf("Comp %d value %f\n",ix,vv);
975 if (vv > 0.0001) {
976 oerr += tot + we * vv;
977 printf("Added %f + %f to give oerr %f\n",tot,we * vv,oerr);
978 }
979 tot += we;
980 }
981 oerr /= tot;
982 if (s->pass == 0)
983 oerr *= 2000.0;
984 else
985 oerr *= 1.0;
986 printf("Final after div by %f oerr = %f\n",tot,oerr);
987
988 printf("~1 Order error %f\n",oerr);
989 rv += oerr;
990 }
991
992 printf("~1 Returning total error %f\n",rv);
993 return rv;
994 }
995
996
997 /* Do a reverse lookup on the mpp */
mpp_rev(mpp * mppo,double limit,double * out,double * in)998 static void mpp_rev(
999 mpp *mppo,
1000 double limit, /* Ink limit */
1001 double *out, /* Device value */
1002 double *in /* Lab target */
1003 ) {
1004 int i, j;
1005 inkmask imask; /* Device Ink mask */
1006 int inn;
1007 revlus rs; /* Reverse lookup structure */
1008 double sr[MAX_CHAN]; /* Search radius */
1009 double tt;
1010 /* !!! This needs to be cached elsewhere !!!! */
1011 static saent *start = NULL; /* Start array */
1012 static int nisay = 0; /* Number in start array */
1013
1014 mppo->get_info(mppo, &imask, &inn, NULL, NULL, NULL, NULL, NULL, NULL);
1015
1016 rs.di = inn; /* Number of device channels */
1017
1018 rs.Lab[0] = in[0]; /* Target PCS value */
1019 rs.Lab[1] = in[1];
1020 rs.Lab[2] = in[2];
1021
1022 rs.dev2lab = mppo->lookup; /* Dev->PCS Lookup function and context */
1023 rs.d2lcntx = (void *)mppo;
1024
1025 rs.ilimit = limit; /* Total ink limit */
1026
1027 {
1028 double Labw[3]; /* Lab value of white */
1029 double Lab[MAX_CHAN][3]; /* Lab value of device primaries */
1030 double min, max;
1031
1032 /* Lookup the L value of all the device primaries */
1033 for (j = 0; j < inn; j++)
1034 out[j] = 0.0;
1035
1036 mppo->lookup(mppo, Labw, out);
1037
1038 for (i = 0; i < inn; i++) {
1039 double tt;
1040 double de;
1041
1042 out[i] = 1.0;
1043 mppo->lookup(mppo, Lab[i], out);
1044
1045 /* Use DE measure heavily weighted towards L only */
1046 tt = L_WEIGHT * (Labw[0] - Lab[i][0]);
1047 de = tt * tt;
1048 tt = 0.4 * (Labw[1] - Lab[i][1]);
1049 de += tt * tt;
1050 tt = 0.2 * (Labw[2] - Lab[i][2]);
1051 de += tt * tt;
1052 rs.oweight[i] = sqrt(de);
1053 out[i] = 0.0;
1054 }
1055
1056 /* Normalise weights from 0 .. 1.0 */
1057 min = 1e6, max = 0.0;
1058 for (j = 0; j < inn; j++) {
1059 if (rs.oweight[j] < min)
1060 min = rs.oweight[j];
1061 if (rs.oweight[j] > max)
1062 max = rs.oweight[j];
1063 }
1064 for (j = 0; j < inn; j++)
1065 rs.oweight[j] = (rs.oweight[j] - min)/(max - min);
1066
1067 {
1068 for (j = 0; j < inn; j++)
1069 rs.sord[j] = j;
1070
1071 for (i = 0; i < (inn-1); i++) {
1072 for (j = i+1; j < inn; j++) {
1073 if (rs.oweight[rs.sord[i]] > rs.oweight[rs.sord[j]]) {
1074 int xx;
1075 xx = rs.sord[i];
1076 rs.sord[i] = rs.sord[j];
1077 rs.sord[j] = xx;
1078 }
1079 }
1080 }
1081 }
1082
1083 for (j = 0; j < inn; j++)
1084 printf("~1 oweight[%d] = %f\n",j,rs.oweight[j]);
1085 for (j = 0; j < inn; j++)
1086 printf("~1 sorted oweight[%d] = %f\n",j,rs.oweight[rs.sord[j]]);
1087 }
1088
1089 /* Initialise the start point array */
1090 if (start == NULL) {
1091 int mxstart;
1092 int steps = 4;
1093
1094 DCOUNT(dix, MAX_CHAN, inn, 0, 0, steps);
1095
1096 printf("~1 initing start point array\n");
1097 for (mxstart = 1, j = 0; j < inn; j++) /* Compute maximum entries */
1098 mxstart *= steps;
1099
1100 printf("~1 mxstart = %d\n",mxstart);
1101 if ((start = malloc(sizeof(saent) * mxstart)) == NULL)
1102 error("mpp_rev malloc of start array failed\n");
1103
1104 nisay = 0;
1105 DC_INIT(dix);
1106
1107 while(!DC_DONE(dix)) {
1108 double sum = 0.0;
1109
1110 /* Figure device values */
1111 for (j = 0; j < inn; j++) {
1112 sum += start[nisay].dev[j] = dix[j]/(steps-1.0);
1113 }
1114
1115 if (sum <= limit) { /* Within ink limit */
1116 double oerr;
1117 double tot;
1118
1119 /* Compute Lab */
1120 mppo->lookup(mppo, start[nisay].Lab, start[nisay].dev);
1121
1122 /* Compute order error */
1123 /* Skip first 3 colorants */
1124 oerr = tot = 0.0;
1125 for (j = 3; j < inn; j++) {
1126 int ix = rs.sord[j]; /* Sorted order index */
1127 double vv = start[nisay].dev[ix];
1128 double we = (double)j - 2.0; /* Increasing weight */
1129
1130 if (vv > 0.0001) {
1131 oerr += tot + we * vv;
1132 }
1133 tot += we;
1134 }
1135 oerr /= tot;
1136 start[nisay].oerr = oerr;
1137
1138 nisay++;
1139 }
1140
1141 DC_INC(dix);
1142 }
1143 printf("~1 start point array done, %d out of %d valid\n",nisay,mxstart);
1144 }
1145
1146 /* Search the start array for closest matching point */
1147 {
1148 double bde = 1e38;
1149 int bix = 0;
1150
1151 for (i = 0; i < nisay; i++) {
1152 double de;
1153
1154 /* Compute delta E */
1155 for (de = 0.0, j = 0; j < 3; j++) {
1156 double tt = rs.Lab[j] - start[i].Lab[j];
1157 de += tt * tt;
1158 }
1159 de += 0.0 * start[i].oerr;
1160 if (de < bde) {
1161 bde = de;
1162 bix = i;
1163 }
1164 }
1165
1166 printf("Start point at index %d, bde = %f, dev = ",bix,bde);
1167 for (j = 0; j < inn; j++) {
1168 printf("%f",start[bix].dev[j]);
1169 if (j < (inn-1))
1170 printf(" ");
1171 }
1172 printf("\n");
1173
1174 for (j = 0; j < inn; j++) {
1175 out[j] = start[bix].dev[j];
1176 sr[j] = 0.1;
1177 }
1178 }
1179
1180 #ifdef NEVER
1181 out[0] = 0.0;
1182 out[1] = 0.0;
1183 out[2] = 0.45;
1184 out[3] = 0.0;
1185 out[4] = 0.0;
1186 out[5] = 0.0;
1187 out[6] = 0.6;
1188 out[7] = 1.0;
1189 #endif
1190
1191 #ifdef NEVER
1192 out[0] = 1.0;
1193 out[1] = 0.0;
1194 out[2] = 0.0;
1195 out[3] = 0.0;
1196 out[4] = 0.0;
1197 out[5] = 0.0;
1198 out[6] = 0.0;
1199 out[7] = 0.0;
1200 #endif
1201
1202 #ifdef NEVER
1203 rs.pass = 0;
1204 if (powell(&tt, inn, out, sr, 0.001, 5000, revoptfunc, (void *)&rs) != 0) {
1205 error("Powell failed inside mpp_rev()");
1206 }
1207 printf("\n\n\n\n\n\n#############################################\n");
1208 printf("~1 after first pass got ");
1209 for (j = 0; j < inn; j++) {
1210 printf("%f",out[j]);
1211 if (j < (inn-1))
1212 printf(" ");
1213 }
1214 printf("\n");
1215 printf("#############################################\n\n\n\n\n\n\n\n");
1216 #endif
1217
1218 #ifndef NEVER
1219 out[0] = 0.0;
1220 out[1] = 0.0;
1221 out[2] = 0.0;
1222 out[3] = 0.0;
1223 out[4] = 0.0;
1224 out[5] = 1.0;
1225 out[6] = 0.0;
1226 out[7] = 0.0;
1227 #endif
1228 #ifndef NEVER
1229 rs.pass = 1;
1230 if (powell(&tt, inn, out, sr, 0.00001, 5000, revoptfunc, (void *)&rs, NULL, NULL) != 0) {
1231 error("Powell failed inside mpp_rev()");
1232 }
1233 #endif
1234
1235 snap2gamut(&rs, out);
1236 }
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255