summaryrefslogtreecommitdiff
path: root/xicc/mpplu.c
diff options
context:
space:
mode:
Diffstat (limited to 'xicc/mpplu.c')
-rw-r--r--xicc/mpplu.c1355
1 files changed, 1355 insertions, 0 deletions
diff --git a/xicc/mpplu.c b/xicc/mpplu.c
new file mode 100644
index 0000000..dcbd1c3
--- /dev/null
+++ b/xicc/mpplu.c
@@ -0,0 +1,1355 @@
+
+/*
+ * Model Printer Profile Lookup test utility.
+ *
+ * Author: Graeme W. Gill
+ * Date: 2002/12/30
+ * Version: 1.00
+ *
+ * Copyright 2002, 2003 Graeme W. Gill
+ * All rights reserved.
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ *
+ */
+
+/* TTBD:
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <fcntl.h>
+#include <string.h>
+#include <math.h>
+#include "numlib.h"
+#include "xicc.h"
+#include "counters.h"
+
+void usage(void) {
+ fprintf(stderr,"Translate colors through an MPP profile, V1.00\n");
+ fprintf(stderr,"Author: Graeme W. Gill\n");
+ fprintf(stderr,"usage: mpplu [-v] [-f func] [-i intent] [-o order] profile\n");
+ fprintf(stderr," -v Verbose\n");
+ fprintf(stderr," -f function f = forward, b = backwards\n");
+ fprintf(stderr," -p oride x = XYZ_PCS, l = Lab_PCS, y = Yxy, s = spectral,\n");
+ fprintf(stderr," -l limit override default ink limit, 1 - N00%%\n");
+ fprintf(stderr," -i illum Choose illuminant for print/transparency spectral data:\n");
+ fprintf(stderr," A, C, D50 (def.), D50M2, D65, F5, F8, F10 or file.sp\n");
+ fprintf(stderr," -o observ Choose CIE Observer for spectral data:\n");
+ fprintf(stderr," 1931_2 (def), 1964_10, S&B 1955_2, shaw, J&V 1978_2\n");
+ fprintf(stderr," -u Use Fluorescent Whitening Agent compensation\n");
+ fprintf(stderr," -g Create gamut output\n");
+ fprintf(stderr," -w Create gamut VRML as well\n");
+ fprintf(stderr," -n Don't add VRML axes\n");
+ fprintf(stderr," -a n Gamut transparency level\n");
+ fprintf(stderr," -d n Gamut surface detail level\n");
+ fprintf(stderr," -t num Invoke debugging test code \"num\" 1..n\n");
+ fprintf(stderr," 1 - check partial derivative for device input\n");
+ fprintf(stderr," 2 - create overlap diagnostic VRML gamut surface\n");
+ fprintf(stderr,"\n");
+ fprintf(stderr," The colors to be translated should be fed into stdin,\n");
+ fprintf(stderr," one input color per line, white space separated.\n");
+ fprintf(stderr," A line starting with a # will be ignored.\n");
+ fprintf(stderr," A line not starting with a number will terminate the program.\n");
+ exit(1);
+}
+
+static void diag_gamut(mpp *p, double gamres, int doaxes, double trans, char *outname);
+static void mpp_rev(mpp *mppo, double limit, double *out, double *in);
+
+int
+main(int argc, char *argv[]) {
+ int fa,nfa; /* argument we're looking at */
+ char prof_name[100];
+ mpp *mppo;
+ int verb = 0;
+ int test = 0; /* special test code */
+ int dogam = 0; /* Create gamut */
+ int dowrl = 0; /* Create VRML gamut */
+ int doaxes = 1; /* Create VRML axes */
+ double trans = 0.0; /* Transparency */
+ double gamres = 0.0; /* Gamut resolution */
+ int repYxy = 0; /* Report Yxy */
+ int repSpec = 0; /* Report Spectral */
+ int bwd = 0; /* Do reverse lookup */
+ double dlimit; /* Device ink limit */
+ double limit = -1.0; /* Used ink limit */
+ int display = 0; /* NZ if display type */
+ int spec = 0; /* Use spectral data flag */
+ int spec_n; /* Number of spectral bands, 0 if not valid */
+ double spec_wl_short; /* First reading wavelength in nm (shortest) */
+ double spec_wl_long; /* Last reading wavelength in nm (longest) */
+ int fwacomp = 0; /* FWA compensation */
+ icxIllumeType illum = icxIT_default; /* Spectral defaults */
+ xspect cust_illum; /* Custom illumination spectrum */
+ icxObserverType observ = icxOT_default;
+ char buf[200];
+ double in[MAX_CHAN], out[MAX_CHAN];
+ int rv = 0;
+
+ inkmask imask; /* Device Ink mask */
+ char *ident = NULL; /* Device colorspec description */
+ icColorSpaceSignature pcss; /* Type of PCS space */
+ int devn, pcsn; /* Number of components */
+
+ icColorSpaceSignature pcsor = icSigLabData; /* Default */
+
+ error_program = argv[0];
+
+ if (argc < 2)
+ usage();
+
+ /* Process the arguments */
+ for(fa = 1;fa < argc;fa++) {
+ nfa = fa; /* skip to nfa if next argument is used */
+ if (argv[fa][0] == '-') { /* Look for any flags */
+ char *na = NULL; /* next argument after flag, null if none */
+
+ if (argv[fa][2] != '\000')
+ na = &argv[fa][2]; /* next is directly after flag */
+ else {
+ if ((fa+1) < argc) {
+ if (argv[fa+1][0] != '-') {
+ nfa = fa + 1;
+ na = argv[nfa]; /* next is seperate non-flag argument */
+ }
+ }
+ }
+
+ if (argv[fa][1] == '?')
+ usage();
+
+ /* Verbosity */
+ else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') {
+ verb = 1;
+ }
+ /* function */
+ else if (argv[fa][1] == 'f' || argv[fa][1] == 'F') {
+ fa = nfa;
+ if (na == NULL) usage();
+ switch (na[0]) {
+ case 'f':
+ case 'F':
+ bwd = 0;
+ break;
+ case 'b':
+ case 'B':
+ bwd = 1;
+ break;
+ default:
+ usage();
+ }
+ }
+
+ /* PCS override */
+ else if (argv[fa][1] == 'p' || argv[fa][1] == 'P') {
+ fa = nfa;
+ if (na == NULL) usage();
+ switch (na[0]) {
+ case 'x':
+ case 'X':
+ pcsor = icSigXYZData;
+ repYxy = repSpec = 0;
+ break;
+ case 'l':
+ case 'L':
+ pcsor = icSigLabData;
+ repYxy = repSpec = 0;
+ break;
+ case 'y':
+ case 'Y':
+ pcsor = icSigXYZData;
+ repYxy = 1;
+ repSpec = 0;
+ break;
+ case 's':
+ case 'S':
+ pcsor = icSigXYZData;
+ repYxy = 0;
+ repSpec = 1;
+ spec = 1;
+ break;
+ default:
+ usage();
+ }
+ }
+
+ /* Ink Limit */
+ else if (argv[fa][1] == 'l' || argv[fa][1] == 'L') {
+ fa = nfa;
+ if (na == NULL) usage();
+ limit = atof(na);
+ }
+
+ /* Spectral Illuminant type */
+ else if (argv[fa][1] == 'i' || argv[fa][1] == 'I') {
+ fa = nfa;
+ if (na == NULL) usage();
+ if (strcmp(na, "A") == 0) {
+ spec = 1;
+ illum = icxIT_A;
+ } else if (strcmp(na, "C") == 0) {
+ spec = 1;
+ illum = icxIT_C;
+ } else if (strcmp(na, "D50") == 0) {
+ spec = 1;
+ illum = icxIT_D50;
+ } else if (strcmp(na, "D50M2") == 0) {
+ spec = 1;
+ illum = icxIT_D50M2;
+ } else if (strcmp(na, "D65") == 0) {
+ spec = 1;
+ illum = icxIT_D65;
+ } else if (strcmp(na, "F5") == 0) {
+ spec = 1;
+ illum = icxIT_F5;
+ } else if (strcmp(na, "F8") == 0) {
+ spec = 1;
+ illum = icxIT_F8;
+ } else if (strcmp(na, "F10") == 0) {
+ spec = 1;
+ illum = icxIT_F10;
+ } else { /* Assume it's a filename */
+ spec = 1;
+ illum = icxIT_custom;
+ if (read_xspect(&cust_illum, na) != 0)
+ usage();
+ }
+ }
+
+ /* Spectral Observer type */
+ else if (argv[fa][1] == 'o' || argv[fa][1] == 'O') {
+ fa = nfa;
+ if (na == NULL) usage();
+ if (strcmp(na, "1931_2") == 0) { /* Classic 2 degree */
+ spec = 1;
+ observ = icxOT_CIE_1931_2;
+ } else if (strcmp(na, "1964_10") == 0) { /* Classic 10 degree */
+ spec = 1;
+ observ = icxOT_CIE_1964_10;
+ } else if (strcmp(na, "1955_2") == 0) { /* Stiles and Burch 1955 2 degree */
+ spec = 1;
+ observ = icxOT_Stiles_Burch_2;
+ } else if (strcmp(na, "1978_2") == 0) { /* Judd and Voss 1978 2 degree */
+ spec = 1;
+ observ = icxOT_Judd_Voss_2;
+ } else if (strcmp(na, "shaw") == 0) { /* Shaw and Fairchilds 1997 2 degree */
+ spec = 1;
+ observ = icxOT_Shaw_Fairchild_2;
+ } else
+ usage();
+ }
+
+ /* Fluerescent Whitner compensation */
+ else if (argv[fa][1] == 'u' || argv[fa][1] == 'U')
+ fwacomp = 1;
+
+ /* Gamut plot */
+ else if (argv[fa][1] == 'g' || argv[fa][1] == 'G')
+ dogam = 1;
+
+ /* VRML plot */
+ else if (argv[fa][1] == 'w' || argv[fa][1] == 'W') {
+ dogam = 1;
+ dowrl = 1;
+ }
+
+ /* No VRML axes */
+ else if (argv[fa][1] == 'n' || argv[fa][1] == 'N') {
+ doaxes = 0;
+ }
+
+ /* Transparency */
+ else if (argv[fa][1] == 'a' || argv[fa][1] == 'A') {
+ fa = nfa;
+ if (na == NULL) usage();
+ trans = atof(na);
+ }
+
+ /* Surface Detail */
+ else if (argv[fa][1] == 'd' || argv[fa][1] == 'D') {
+ fa = nfa;
+ if (na == NULL) usage();
+ gamres = atof(na);
+ dogam = 1;
+ }
+
+ /* Test code */
+ else if (argv[fa][1] == 't' || argv[fa][1] == 'T') {
+ fa = nfa;
+ if (na == NULL) usage();
+ test = atoi(na);
+ }
+
+ else
+ usage();
+ } else
+ break;
+ }
+
+ if (fa >= argc || argv[fa][0] == '-') usage();
+ strcpy(prof_name,argv[fa]);
+
+ if ((mppo = new_mpp()) == NULL)
+ error ("Creation of MPP object failed");
+
+ if ((rv = mppo->read_mpp(mppo,prof_name)) != 0)
+ error ("%d, %s",rv,mppo->err);
+
+ mppo->get_info(mppo, &imask, &devn, &dlimit, &spec_n, &spec_wl_short, &spec_wl_long, NULL, &display);
+ ident = icx_inkmask2char(imask, 1);
+
+ if (verb) {
+ printf("MPP profile with %d colorants, type %s, TAC %f\n",devn,ident, dlimit);
+ if (display)
+ printf("MPP profile is for a display type device\n");
+ }
+
+ if (limit <= 0.0 || dlimit < limit)
+ limit = dlimit;
+
+ pcss = pcsor;
+ pcsn = 3;
+
+ if (spec && spec_n == 0) {
+ error("Spectral profile needed for spectral result, custom illuminant, observer or FWA");
+ }
+
+ /* Select CIE return value details */
+ if ((rv = mppo->set_ilob(mppo, illum, &cust_illum, observ, NULL, pcss, fwacomp)) != 0) {
+ if (rv == 1)
+ error("Spectral profile needed for custom illuminant, observer or FWA");
+ error("Error setting illuminant, observer, or FWA");
+ }
+
+ if (test != 0) {
+ printf("!!!!! Running special test code no %d !!!!!\n",test);
+
+ if (test == 1) {
+ double **dv, **rdv;
+
+ dv = dmatrix(0, pcsn-1, 0, devn-1);
+ rdv = dmatrix(0, pcsn-1, 0, devn-1);
+
+ printf("Checking partial derivative at each input value\n");
+ /* Process colors to translate */
+ for (;;) {
+ int i,j;
+ char *bp, *nbp;
+ double tout[MAX_CHAN];
+
+ /* Read in the next line */
+ if (fgets(buf, 200, stdin) == NULL)
+ break;
+ if (buf[0] == '#') {
+ fprintf(stdout,"%s\n",buf);
+ continue;
+ }
+ /* For each input number */
+ for (bp = buf-1, nbp = buf, i = 0; i < MAX_CHAN; i++) {
+ bp = nbp;
+ in[i] = strtod(bp, &nbp);
+ if (nbp == bp)
+ break; /* Failed */
+ }
+ if (i == 0)
+ break;
+
+ /* Do conversion */
+ mppo->lookup(mppo, out, in);
+ mppo->dlookup(mppo, out, dv, in);
+
+ for (j = 0; j < devn; j++) {
+ double del = 1e-9;
+
+ if (in[j] > 0.5)
+ del = -del;
+
+ in[j] += del;
+ mppo->lookup(mppo, tout, in);
+ in[j] -= del;
+
+ for (i = 0; i < pcsn; i++) {
+ rdv[i][j] = (tout[i] - out[i])/del;
+ }
+ }
+
+ /* Output the results */
+ for (j = 0; j < devn; j++) {
+ if (j > 0)
+ fprintf(stdout," %f",in[j]);
+ else
+ fprintf(stdout,"%f",in[j]);
+ }
+ printf(" [%s] -> ", ident);
+
+ for (j = 0; j < pcsn; j++) {
+ if (j > 0)
+ fprintf(stdout," %f",out[j]);
+ else
+ fprintf(stdout,"%f",out[j]);
+ }
+ printf(" [%s]\n", icm2str(icmColorSpaceSignature, pcss));
+
+ /* Print the derivatives */
+ for (i = 0; i < pcsn; i++) {
+
+ printf("Output chan %d: ",i);
+ for (j = 0; j < devn; j++) {
+ if (j < (devn-1))
+ fprintf(stdout,"%f ref %f, ",dv[i][j], rdv[i][j]);
+ else
+ fprintf(stdout,"%f ref %f\n",dv[i][j], rdv[i][j]);
+ }
+ }
+ }
+
+ free_dmatrix(dv, 0, pcsn-1, 0, devn-1);
+ free_dmatrix(rdv, 0, pcsn-1, 0, devn-1);
+
+ } else if (test == 2) {
+ char *xl, gam_name[100];
+
+ strcpy(gam_name, prof_name);
+ if ((xl = strrchr(gam_name, '.')) == NULL) /* Figure where extention is */
+ xl = gam_name + strlen(gam_name);
+
+ strcpy(xl,".wrl");
+ diag_gamut(mppo, gamres, doaxes, trans, gam_name);
+
+ } else {
+ printf("Unknown test!\n");
+ }
+
+ } else if (dogam) {
+ gamut *gam;
+ char *xl, gam_name[100];
+ int docusps = 1;
+
+ if ((gam = mppo->get_gamut(mppo, gamres)) == NULL)
+ error("get_gamut failed\n");
+
+ strcpy(gam_name, prof_name);
+ if ((xl = strrchr(gam_name, '.')) == NULL) /* Figure where extention is */
+ xl = gam_name + strlen(gam_name);
+
+ strcpy(xl,".gam");
+ if (gam->write_gam(gam,gam_name))
+ error ("write gamut failed on '%s'",gam_name);
+
+ if (dowrl) {
+ strcpy(xl,".wrl");
+ if (gam->write_vrml(gam,gam_name, doaxes, docusps))
+ error ("write vrml failed on '%s'",gam_name);
+ }
+
+ gam->del(gam);
+
+ } else {
+ /* Normal color lookup */
+
+ if (repYxy) { /* report Yxy rather than XYZ */
+ if (pcss == icSigXYZData)
+ pcss = icSigYxyData;
+ }
+
+ /* Process colors to translate */
+ for (;;) {
+ int i,j;
+ char *bp, *nbp;
+
+ /* Read in the next line */
+ if (fgets(buf, 200, stdin) == NULL)
+ break;
+ if (buf[0] == '#') {
+ fprintf(stdout,"%s\n",buf);
+ continue;
+ }
+ /* For each input number */
+ for (bp = buf-1, nbp = buf, i = 0; i < MAX_CHAN; i++) {
+ bp = nbp;
+ in[i] = strtod(bp, &nbp);
+ if (nbp == bp)
+ break; /* Failed */
+ }
+ if (i == 0)
+ break;
+
+ if (!bwd) {
+
+ if (repSpec) {
+ xspect ospec;
+
+ /* Do lookup of spectrum */
+ mppo->lookup_spec(mppo, &ospec, in);
+
+ /* Output the results */
+ for (j = 0; j < devn; j++) {
+ if (j > 0)
+ fprintf(stdout," %f",in[j]);
+ else
+ fprintf(stdout,"%f",in[j]);
+ }
+ printf(" [%s] -> ", ident);
+
+ for (j = 0; j < spec_n; j++) {
+ if (j > 0)
+ fprintf(stdout," %f",ospec.spec[j]);
+ else
+ fprintf(stdout,"%f",ospec.spec[j]);
+ }
+
+ printf(" [%3.0f .. %3.0f nm]\n", spec_wl_short, spec_wl_long);
+
+ } else {
+
+ /* Do conversion */
+ mppo->lookup(mppo, out, in);
+
+ /* Output the results */
+ for (j = 0; j < devn; j++) {
+ if (j > 0)
+ fprintf(stdout," %f",in[j]);
+ else
+ fprintf(stdout,"%f",in[j]);
+ }
+ printf(" [%s] -> ", ident);
+
+ if (repYxy && pcss == icSigYxyData) {
+ double X = out[0];
+ double Y = out[1];
+ double Z = out[2];
+ double sum = X + Y + Z;
+ if (sum < 1e-6) {
+ out[0] = out[1] = out[2] = 0.0;
+ } else {
+ out[0] = Y;
+ out[1] = X/sum;
+ out[2] = Y/sum;
+ }
+ }
+ for (j = 0; j < pcsn; j++) {
+ if (j > 0)
+ fprintf(stdout," %f",out[j]);
+ else
+ fprintf(stdout,"%f",out[j]);
+ }
+
+ printf(" [%s]\n", icm2str(icmColorSpaceSignature, pcss));
+ }
+
+ } else { /* Do a reverse lookup */
+
+ if (repYxy && pcss == icSigYxyData) {
+ double Y = in[0];
+ double x = in[1];
+ double y = in[2];
+ double z = 1.0 - x - y;
+ double sum;
+ if (y < 1e-6) {
+ in[0] = in[1] = in[2] = 0.0;
+ } else {
+ sum = Y/y;
+ in[0] = x * sum;
+ in[1] = Y;
+ in[2] = z * sum;
+ }
+ }
+
+ /* Do conversion */
+ mpp_rev(mppo, limit, out, in);
+
+ /* Output the results */
+ for (j = 0; j < pcsn; j++) {
+ if (j > 0)
+ fprintf(stdout," %f",in[j]);
+ else
+ fprintf(stdout,"%f",in[j]);
+ }
+ printf(" [%s] -> ", icm2str(icmColorSpaceSignature, pcss));
+
+ for (j = 0; j < devn; j++) {
+ if (j > 0)
+ fprintf(stdout," %f",out[j]);
+ else
+ fprintf(stdout,"%f",out[j]);
+ }
+
+ printf(" [%s]\n", ident);
+ }
+ }
+ }
+
+ free(ident);
+ mppo->del(mppo);
+
+ return 0;
+}
+
+
+/* -------------------------------------------- */
+/* Code for special gamut surface plot */
+
+#define GAMUT_LCENT 50
+
+/* Create a diagnostic gamut, illustrating */
+/* device space "fold-over" */
+/* (This will be in the current PCS, but assumed to be Lab) */
+static void diag_gamut(
+mpp *p, /* This */
+double detail, /* Gamut resolution detail */
+int doaxes,
+double trans, /* Transparency */
+char *outname /* Output VRML file */
+) {
+ int i, j;
+ FILE *wrl;
+ struct {
+ double x, y, z;
+ double wx, wy, wz;
+ double r, g, b;
+ } axes[5] = {
+ { 0, 0, 50-GAMUT_LCENT, 2, 2, 100, .7, .7, .7 }, /* L axis */
+ { 50, 0, 0-GAMUT_LCENT, 100, 2, 2, 1, 0, 0 }, /* +a (red) axis */
+ { 0, -50, 0-GAMUT_LCENT, 2, 100, 2, 0, 0, 1 }, /* -b (blue) axis */
+ { -50, 0, 0-GAMUT_LCENT, 100, 2, 2, 0, 1, 0 }, /* -a (green) axis */
+ { 0, 50, 0-GAMUT_LCENT, 2, 100, 2, 1, 1, 0 }, /* +b (yellow) axis */
+ };
+ int vix; /* Vertex index */
+ DCOUNT(coa, MAX_CHAN, p->n, 0, 0, 2);
+ double col[MPP_MXCCOMB][3]; /* Color asigned to each major vertex */
+ int res;
+
+ /* Asign some colors to the combination nodes */
+ for (i = 0; i < p->nn; i++) {
+ int a, b, c, j;
+ double h;
+
+ j = (i ^ 0x5a5a5a5a) % p->nn;
+ h = (double)j/(p->nn-1);
+
+ /* Make fully saturated with chosen hue */
+ if (h < 1.0/3.0) {
+ a = 0;
+ b = 1;
+ c = 2;
+ } else if (h < 2.0/3.0) {
+ a = 1;
+ b = 2;
+ c = 0;
+ h -= 1.0/3.0;
+ } else {
+ a = 2;
+ b = 0;
+ c = 1;
+ h -= 2.0/3.0;
+ }
+ h *= 3.0;
+
+ col[i][a] = (1.0 - h);
+ col[i][b] = h;
+ col[i][c] = d_rand(0.0, 1.0);
+ }
+
+ if (detail > 0.0)
+ res = (int)(100.0/detail); /* Establish an appropriate sampling density */
+ else
+ res = 4;
+
+ if (res < 2)
+ res = 2;
+
+ if ((wrl = fopen(outname,"w")) == NULL)
+ error("Error opening wrl output file '%s'",outname);
+
+ /* Spit out a VRML 2 Object surface of gamut */
+ fprintf(wrl,"#VRML V2.0 utf8\n");
+ fprintf(wrl,"\n");
+ fprintf(wrl,"# Created by the Argyll CMS\n");
+ fprintf(wrl,"Transform {\n");
+ fprintf(wrl,"children [\n");
+ fprintf(wrl," NavigationInfo {\n");
+ fprintf(wrl," type \"EXAMINE\" # It's an object we examine\n");
+ fprintf(wrl," } # We'll add our own light\n");
+ fprintf(wrl,"\n");
+ fprintf(wrl," DirectionalLight {\n");
+ fprintf(wrl," direction 0 0 -1 # Light illuminating the scene\n");
+ fprintf(wrl," direction 0 -1 0 # Light illuminating the scene\n");
+ fprintf(wrl," }\n");
+ fprintf(wrl,"\n");
+ fprintf(wrl," Viewpoint {\n");
+ fprintf(wrl," position 0 0 340 # Position we view from\n");
+ fprintf(wrl," }\n");
+ fprintf(wrl,"\n");
+ if (doaxes != 0) {
+ fprintf(wrl,"# Lab axes as boxes:\n");
+ for (i = 0; i < 5; i++) {
+ fprintf(wrl,"Transform { translation %f %f %f\n", axes[i].x, axes[i].y, axes[i].z);
+ fprintf(wrl,"\tchildren [\n");
+ fprintf(wrl,"\t\tShape{\n");
+ fprintf(wrl,"\t\t\tgeometry Box { size %f %f %f }\n",
+ axes[i].wx, axes[i].wy, axes[i].wz);
+ fprintf(wrl,"\t\t\tappearance Appearance { material Material ");
+ fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[i].r, axes[i].g, axes[i].b);
+ fprintf(wrl,"\t\t}\n");
+ fprintf(wrl,"\t]\n");
+ fprintf(wrl,"}\n");
+ }
+ fprintf(wrl,"\n");
+ }
+ fprintf(wrl," Transform {\n");
+ fprintf(wrl," translation 0 0 0\n");
+ fprintf(wrl," children [\n");
+ fprintf(wrl," Shape { \n");
+ fprintf(wrl," geometry IndexedFaceSet {\n");
+ fprintf(wrl," solid FALSE\n"); /* Don't back face cull */
+ fprintf(wrl," convex TRUE\n");
+ fprintf(wrl,"\n");
+ fprintf(wrl," coord Coordinate { \n");
+ fprintf(wrl," point [ # Verticy coordinates\n");
+
+
+ /* Itterate over all the faces in the device space */
+ /* generating the vertx positions. */
+ DC_INIT(coa);
+ vix = 0;
+ while(!DC_DONE(coa)) {
+ int e, m1, m2;
+ double in[MAX_CHAN];
+ double out[3];
+ double sum;
+
+ /* Scan only device surface */
+ for (m1 = 0; m1 < p->n; m1++) {
+ if (coa[m1] != 0)
+ continue;
+
+ for (m2 = m1 + 1; m2 < p->n; m2++) {
+ int x, y;
+
+ if (coa[m2] != 0)
+ continue;
+
+ for (sum = 0.0, e = 0; e < p->n; e++)
+ in[e] = (double)coa[e]; /* Base value */
+
+ /* Scan over 2D device space face */
+ for (x = 0; x < res; x++) { /* step over surface */
+ in[m1] = x/(res - 1.0);
+ for (y = 0; y < res; y++) {
+ in[m2] = y/(res - 1.0);
+
+ p->lookup(p, out, in);
+ fprintf(wrl,"%f %f %f,\n",out[1], out[2], out[0]-50.0);
+ vix++;
+ }
+ }
+ }
+ }
+ /* Increment index within block */
+ DC_INC(coa);
+ }
+
+ fprintf(wrl," ]\n");
+ fprintf(wrl," }\n");
+ fprintf(wrl,"\n");
+ fprintf(wrl," coordIndex [ # Indexes of poligon Verticies \n");
+
+ /* Itterate over all the faces in the device space */
+ /* generating the quadrilateral indexes. */
+ DC_INIT(coa);
+ vix = 0;
+ while(!DC_DONE(coa)) {
+ int e, m1, m2;
+ double in[MAX_CHAN];
+ double sum;
+
+ /* Scan only device surface */
+ for (m1 = 0; m1 < p->n; m1++) {
+ if (coa[m1] != 0)
+ continue;
+
+ for (m2 = m1 + 1; m2 < p->n; m2++) {
+ int x, y;
+
+ if (coa[m2] != 0)
+ continue;
+
+ for (sum = 0.0, e = 0; e < p->n; e++)
+ in[e] = (double)coa[e]; /* Base value */
+
+ /* Scan over 2D device space face */
+ for (x = 0; x < res; x++) { /* step over surface */
+ for (y = 0; y < res; y++) {
+
+ if (x < (res-1) && y < (res-1)) {
+ fprintf(wrl,"%d, %d, %d, %d, -1\n",
+ vix, vix + 1, vix + 1 + res, vix + res);
+ }
+ vix++;
+ }
+ }
+ }
+ }
+ /* Increment index within block */
+ DC_INC(coa);
+ }
+
+ fprintf(wrl," ]\n");
+ fprintf(wrl,"\n");
+ fprintf(wrl," colorPerVertex TRUE\n");
+ fprintf(wrl," color Color {\n");
+ fprintf(wrl," color [ # RGB colors of each vertex\n");
+
+ /* Itterate over all the faces in the device space */
+ /* generating the vertx colors. */
+ DC_INIT(coa);
+ vix = 0;
+ while(!DC_DONE(coa)) {
+ int e, m1, m2;
+ double in[MAX_CHAN];
+ double sum;
+
+ /* Scan only device surface */
+ for (m1 = 0; m1 < p->n; m1++) {
+ if (coa[m1] != 0)
+ continue;
+
+ for (m2 = m1 + 1; m2 < p->n; m2++) {
+ int x, y;
+
+ if (coa[m2] != 0)
+ continue;
+
+ for (sum = 0.0, e = 0; e < p->n; e++)
+ in[e] = (double)coa[e]; /* Base value */
+
+ /* Scan over 2D device space face */
+ for (x = 0; x < res; x++) { /* step over surface */
+ double xb = x/(res - 1.0);
+ for (y = 0; y < res; y++) {
+ int v0, v1, v2, v3;
+ double yb = y/(res - 1.0);
+ double rgb[3];
+
+ for (v0 = 0, e = 0; e < p->n; e++)
+ v0 |= coa[e] ? (1 << e) : 0; /* Binary index */
+
+ v1 = v0 | (1 << m2); /* Y offset */
+ v2 = v0 | (1 << m2) | (1 << m1); /* X+Y offset */
+ v3 = v0 | (1 << m1); /* Y offset */
+
+ /* Linear interp between the main verticies */
+ for (j = 0; j < 3; j++) {
+ rgb[j] = (1.0 - yb) * (1.0 - xb) * col[v0][j]
+ + yb * (1.0 - xb) * col[v1][j]
+ + (1.0 - yb) * xb * col[v3][j]
+ + yb * xb * col[v2][j];
+ }
+ fprintf(wrl,"%f %f %f,\n",rgb[1], rgb[2], rgb[0]);
+ vix++;
+ }
+ }
+ }
+ }
+ /* Increment index within block */
+ DC_INC(coa);
+ }
+
+ fprintf(wrl," ] \n");
+ fprintf(wrl," }\n");
+ fprintf(wrl," }\n");
+ fprintf(wrl," appearance Appearance { \n");
+ fprintf(wrl," material Material {\n");
+ fprintf(wrl," transparency %f\n",trans);
+ fprintf(wrl," ambientIntensity 0.3\n");
+ fprintf(wrl," shininess 0.5\n");
+ fprintf(wrl," }\n");
+ fprintf(wrl," }\n");
+ fprintf(wrl," } # end Shape\n");
+ fprintf(wrl," ]\n");
+ fprintf(wrl," }\n");
+
+ fprintf(wrl,"\n");
+ fprintf(wrl," ] # end of children for world\n");
+ fprintf(wrl,"}\n");
+
+ if (fclose(wrl) != 0)
+ error("Error closing output file '%s'",outname);
+}
+
+/* -------------------------------------------- */
+/* Reverse lookup support */
+/* This is for developing the appropriate reverse lookup */
+/* code for xsep.c */
+
+/*
+ * TTBD:
+ * Not handlink rule or separate black
+ * Not handling linearisation callback for ink limit.
+ * Not allowing for other possible secondary limits/goals/tunings.
+ */
+
+/*
+ * Descriptionr:
+ *
+ * The primary reverse lookup optimisation goals are to remain within
+ * the device gamut, and to match the target PCS.
+ *
+ * The secondary optimisation goals for solving this under constrained
+ * problem can be chosen to a achieve a wide variety of possible aims.
+ *
+ * 1) One that is applicable to screened devices might be to use the extra
+ * inks to minimise the visiblility of screening. For screening resolutions
+ * above 50 lpi (equivalent to 100 dpi stocastic screens) only luminance
+ * contrast will be relavant, so priority to the inks closest to paper white
+ * measured purely by L distance is appropriate. (In practice I've used
+ * a measure that adds a small color distance component.)
+ *
+ * 2) Another possible goal would be to optimise fit to a desired spectral
+ * profile, if spectral information is avaiable as an aim. The goal would
+ * be best fit weighted to the visual sensitivity curve.
+ *
+ * 3) Another possible secondary goal might be to minimise the total
+ * amount of ink used, to minimise cost, drying time, media ink loading.
+ *
+ * 4) A fourth possible secondary goal might be to choose ink combinations
+ * that are the most robust given an error in device values.
+ *
+ * In practice these secondary goals need not be entirely exclusive,
+ * but the overall secondary goal could be a weighted blending between
+ * these goals. Overall the ink limit (TAC) is extremely important,
+ * since this will be the primary thing that stops large amounds of
+ * light ink being used at all times.
+ *
+ * Numerous other tweaks, limits or goals (such as secondary combination
+ * ink limits, exclusions such as Cyan/Oraange, Magenta/Green)
+ * could also be applied in the reverse optimisation routine.
+ *
+ */
+
+
+#ifdef NEVER
+/* These weights are the "theoreticaly correct" weightings, since */
+/* at 50 lpi or higher, the color contrast sensitivity should be close to 0 */
+# define L_WEIGHT 1.0
+# define a_WEIGHT 0.0
+# define b_WEIGHT 0.0
+#else
+/* These weights give us our "expected" ink ordering of */
+/* Yellow, light Cyan/Magenta, Orange/Green, Cyan/Magenta, Black. */
+# define L_WEIGHT 1.0
+# define a_WEIGHT 0.4
+# define b_WEIGHT 0.2
+#endif
+
+/* Start array entry */
+typedef struct {
+ double dev[MAX_CHAN]; /* Device value */
+ double Lab[3]; /* Lab value */
+ double oerr; /* Order error */
+} saent;
+
+/* Context for reverse lookup */
+typedef struct {
+ int pass;
+ int di; /* Number of device channels */
+ double Lab[3]; /* Lab target value */
+ void (*dev2lab) (mpp *d2lcntx, double *out, double *in); /* Device to Lab callback */
+ mpp *d2lcntx; /* dev2lab callback context */
+ double ilimit; /* Total limit */
+ int sord[MAX_CHAN]; /* Sorted order index */
+ double oweight[MAX_CHAN]; /* Order weighting (not used ?) */
+} revlus;
+
+/* Return the largest distance of the point outside the device gamut. */
+/* This will be 0 if inside the gamut, and > 0 if outside. */
+static double
+dist2gamut(revlus *s, double *d) {
+ int e;
+ int di = s->di;
+ double tt, dd = 0.0;
+ double ss = 0.0;
+
+ for (e = 0; e < di; e++) {
+ ss += d[e];
+
+ tt = 0.0 - d[e];
+ if (tt > 0.0) {
+ if (tt > dd)
+ dd = tt;
+ }
+ tt = d[e] - 1.0;
+ if (tt > 0.0) {
+ if (tt > dd)
+ dd = tt;
+ }
+ }
+ tt = ss - s->ilimit;
+ if (tt > 0.0) {
+ if (tt > dd)
+ dd = tt;
+ }
+ return dd;
+}
+
+/* Snap a point to the device gamut boundary. */
+/* Return nz if it has been snapped. */
+static int snap2gamut(revlus *s, double *d) {
+ int e;
+ int di = s->di;
+ double dd; /* Smallest distance */
+ double ss; /* Sum */
+ int rv = 0;
+
+ /* Snap to ink limit first */
+ for (ss = 0.0, e = 0; e < di; e++)
+ ss += d[e];
+ dd = fabs(ss - s->ilimit);
+
+ if (dd < 0.0) {
+ int j;
+ for (j = 0; j < di; j++)
+ d[j] *= s->ilimit/ss; /* Snap to ink limit */
+ rv = 1;
+ }
+
+ /* Now snap to any other dimension */
+ for (e = 0; e < di; e++) {
+
+ dd = fabs(d[e] - 0.0);
+ if (dd < 0.0) {
+ d[e] = 0.0; /* Snap to orthogonal boundary */
+ rv = 1;
+ }
+ dd = fabs(1.0 - d[e]);
+ if (dd < 0.0) {
+ d[e] = 1.0; /* Snap to orthogonal boundary */
+ rv = 1;
+ }
+ }
+
+ return rv;
+}
+
+/* Reverse optimisation function handed to powell() */
+static double revoptfunc(void *edata, double *v) {
+ revlus *s = (revlus *)edata;
+ double rv;
+
+printf("~1 target %f %f %f\n",s->Lab[0],s->Lab[1],s->Lab[2]);
+
+ if ((rv = (dist2gamut(s, v))) > 0.0) {
+// rv = rv * 1000.0 + 45000.0; /* Discourage being out of gamut */
+ rv = rv * 5e6; /* Discourage being out of gamut */
+
+ }
+printf("~1 out of gamut error = %f\n",rv);
+ {
+ int j;
+ double oerr;
+ double Lab[3];
+ double tot;
+
+ /* Convert device to Lab */
+ s->dev2lab(s->d2lcntx, Lab, v);
+
+ /* Accumulate total delta E squared */
+ for (j = 0; j < 3; j++) {
+ double tt = s->Lab[j] - Lab[j];
+ rv += tt * tt;
+ }
+
+printf("~1 Delta E squared = %f\n",rv);
+
+ /* Skip first 3 colorants */
+ oerr = tot = 0.0;
+printf("oerr = %f\n",oerr);
+ for (j = 3; j < s->di; j++) {
+ int ix = s->sord[j]; /* Sorted order index */
+ double vv = v[ix];
+ double we = (double)j - 2.0; /* Increasing weight */
+
+printf("Comp %d value %f\n",ix,vv);
+ if (vv > 0.0001) {
+ oerr += tot + we * vv;
+printf("Added %f + %f to give oerr %f\n",tot,we * vv,oerr);
+ }
+ tot += we;
+ }
+ oerr /= tot;
+ if (s->pass == 0)
+ oerr *= 2000.0;
+ else
+ oerr *= 1.0;
+printf("Final after div by %f oerr = %f\n",tot,oerr);
+
+printf("~1 Order error %f\n",oerr);
+ rv += oerr;
+ }
+
+printf("~1 Returning total error %f\n",rv);
+ return rv;
+}
+
+
+/* Do a reverse lookup on the mpp */
+static void mpp_rev(
+mpp *mppo,
+double limit, /* Ink limit */
+double *out, /* Device value */
+double *in /* Lab target */
+) {
+ int i, j;
+ inkmask imask; /* Device Ink mask */
+ int inn;
+ revlus rs; /* Reverse lookup structure */
+ double sr[MAX_CHAN]; /* Search radius */
+ double tt;
+ /* !!! This needs to be cached elsewhere !!!! */
+ static saent *start = NULL; /* Start array */
+ static int nisay = 0; /* Number in start array */
+
+ mppo->get_info(mppo, &imask, &inn, NULL, NULL, NULL, NULL, NULL, NULL);
+
+ rs.di = inn; /* Number of device channels */
+
+ rs.Lab[0] = in[0]; /* Target PCS value */
+ rs.Lab[1] = in[1];
+ rs.Lab[2] = in[2];
+
+ rs.dev2lab = mppo->lookup; /* Dev->PCS Lookup function and context */
+ rs.d2lcntx = (void *)mppo;
+
+ rs.ilimit = limit; /* Total ink limit */
+
+ {
+ double Labw[3]; /* Lab value of white */
+ double Lab[MAX_CHAN][3]; /* Lab value of device primaries */
+ double min, max;
+
+ /* Lookup the L value of all the device primaries */
+ for (j = 0; j < inn; j++)
+ out[j] = 0.0;
+
+ mppo->lookup(mppo, Labw, out);
+
+ for (i = 0; i < inn; i++) {
+ double tt;
+ double de;
+
+ out[i] = 1.0;
+ mppo->lookup(mppo, Lab[i], out);
+
+ /* Use DE measure heavily weighted towards L only */
+ tt = L_WEIGHT * (Labw[0] - Lab[i][0]);
+ de = tt * tt;
+ tt = 0.4 * (Labw[1] - Lab[i][1]);
+ de += tt * tt;
+ tt = 0.2 * (Labw[2] - Lab[i][2]);
+ de += tt * tt;
+ rs.oweight[i] = sqrt(de);
+ out[i] = 0.0;
+ }
+
+ /* Normalise weights from 0 .. 1.0 */
+ min = 1e6, max = 0.0;
+ for (j = 0; j < inn; j++) {
+ if (rs.oweight[j] < min)
+ min = rs.oweight[j];
+ if (rs.oweight[j] > max)
+ max = rs.oweight[j];
+ }
+ for (j = 0; j < inn; j++)
+ rs.oweight[j] = (rs.oweight[j] - min)/(max - min);
+
+ {
+ for (j = 0; j < inn; j++)
+ rs.sord[j] = j;
+
+ for (i = 0; i < (inn-1); i++) {
+ for (j = i+1; j < inn; j++) {
+ if (rs.oweight[rs.sord[i]] > rs.oweight[rs.sord[j]]) {
+ int xx;
+ xx = rs.sord[i];
+ rs.sord[i] = rs.sord[j];
+ rs.sord[j] = xx;
+ }
+ }
+ }
+ }
+
+for (j = 0; j < inn; j++)
+printf("~1 oweight[%d] = %f\n",j,rs.oweight[j]);
+for (j = 0; j < inn; j++)
+printf("~1 sorted oweight[%d] = %f\n",j,rs.oweight[rs.sord[j]]);
+ }
+
+ /* Initialise the start point array */
+ if (start == NULL) {
+ int mxstart;
+ int steps = 4;
+
+ DCOUNT(dix, MAX_CHAN, inn, 0, 0, steps);
+
+printf("~1 initing start point array\n");
+ for (mxstart = 1, j = 0; j < inn; j++) /* Compute maximum entries */
+ mxstart *= steps;
+
+printf("~1 mxstart = %d\n",mxstart);
+ if ((start = malloc(sizeof(saent) * mxstart)) == NULL)
+ error("mpp_rev malloc of start array failed\n");
+
+ nisay = 0;
+ DC_INIT(dix);
+
+ while(!DC_DONE(dix)) {
+ double sum = 0.0;
+
+ /* Figure device values */
+ for (j = 0; j < inn; j++) {
+ sum += start[nisay].dev[j] = dix[j]/(steps-1.0);
+ }
+
+ if (sum <= limit) { /* Within ink limit */
+ double oerr;
+ double tot;
+
+ /* Compute Lab */
+ mppo->lookup(mppo, start[nisay].Lab, start[nisay].dev);
+
+ /* Compute order error */
+ /* Skip first 3 colorants */
+ oerr = tot = 0.0;
+ for (j = 3; j < inn; j++) {
+ int ix = rs.sord[j]; /* Sorted order index */
+ double vv = start[nisay].dev[ix];
+ double we = (double)j - 2.0; /* Increasing weight */
+
+ if (vv > 0.0001) {
+ oerr += tot + we * vv;
+ }
+ tot += we;
+ }
+ oerr /= tot;
+ start[nisay].oerr = oerr;
+
+ nisay++;
+ }
+
+ DC_INC(dix);
+ }
+printf("~1 start point array done, %d out of %d valid\n",nisay,mxstart);
+ }
+
+ /* Search the start array for closest matching point */
+ {
+ double bde = 1e38;
+ int bix = 0;
+
+ for (i = 0; i < nisay; i++) {
+ double de;
+
+ /* Compute delta E */
+ for (de = 0.0, j = 0; j < 3; j++) {
+ double tt = rs.Lab[j] - start[i].Lab[j];
+ de += tt * tt;
+ }
+ de += 0.0 * start[i].oerr;
+ if (de < bde) {
+ bde = de;
+ bix = i;
+ }
+ }
+
+printf("Start point at index %d, bde = %f, dev = ",bix,bde);
+for (j = 0; j < inn; j++) {
+printf("%f",start[bix].dev[j]);
+if (j < (inn-1))
+printf(" ");
+}
+printf("\n");
+
+ for (j = 0; j < inn; j++) {
+ out[j] = start[bix].dev[j];
+ sr[j] = 0.1;
+ }
+ }
+
+#ifdef NEVER
+ out[0] = 0.0;
+ out[1] = 0.0;
+ out[2] = 0.45;
+ out[3] = 0.0;
+ out[4] = 0.0;
+ out[5] = 0.0;
+ out[6] = 0.6;
+ out[7] = 1.0;
+#endif
+
+#ifdef NEVER
+ out[0] = 1.0;
+ out[1] = 0.0;
+ out[2] = 0.0;
+ out[3] = 0.0;
+ out[4] = 0.0;
+ out[5] = 0.0;
+ out[6] = 0.0;
+ out[7] = 0.0;
+#endif
+
+#ifdef NEVER
+ rs.pass = 0;
+ if (powell(&tt, inn, out, sr, 0.001, 5000, revoptfunc, (void *)&rs) != 0) {
+ error("Powell failed inside mpp_rev()");
+ }
+printf("\n\n\n\n\n\n#############################################\n");
+printf("~1 after first pass got ");
+for (j = 0; j < inn; j++) {
+printf("%f",out[j]);
+if (j < (inn-1))
+printf(" ");
+}
+printf("\n");
+printf("#############################################\n\n\n\n\n\n\n\n");
+#endif
+
+#ifndef NEVER
+ out[0] = 0.0;
+ out[1] = 0.0;
+ out[2] = 0.0;
+ out[3] = 0.0;
+ out[4] = 0.0;
+ out[5] = 1.0;
+ out[6] = 0.0;
+ out[7] = 0.0;
+#endif
+#ifndef NEVER
+ rs.pass = 1;
+ if (powell(&tt, inn, out, sr, 0.00001, 5000, revoptfunc, (void *)&rs, NULL, NULL) != 0) {
+ error("Powell failed inside mpp_rev()");
+ }
+#endif
+
+ snap2gamut(&rs, out);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+