diff options
Diffstat (limited to 'xicc/mpplu.c')
-rw-r--r-- | xicc/mpplu.c | 1355 |
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); +} + + + + + + + + + + + + + + + + + + |