diff options
Diffstat (limited to 'xicc/mpp.c')
| -rw-r--r-- | xicc/mpp.c | 4446 |
1 files changed, 4446 insertions, 0 deletions
diff --git a/xicc/mpp.c b/xicc/mpp.c new file mode 100644 index 0000000..02b3019 --- /dev/null +++ b/xicc/mpp.c @@ -0,0 +1,4446 @@ + +/* + * Argyll Color Correction System + * Model Printer Profile object. + * + * Author: Graeme W. Gill + * Date: 24/2/2002 + * + * Copyright 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. + */ + +/* This version (based on mpp_x1) has n * 2^(n-1) shape params, */ +/* used for linear interpolation of the shaping correction. */ + +/* + * This program takes in the scattered test chart + * points, and creates a model based forward printer profile + * (Device -> CIE + spectral), based on heuristic extensions to + * the Neugenbauer model. + * It is designed to handle an arbitrary number of colorants, + * and in the future, (optionaly) create an aproximate ink overlap/mixing model + * to allow synthesis of a forward model for a hyperthetical + * similar printing process with additional inks. + * + * The model is as follows: + * + * A per input channel transfer curve, that models + * per channel dot gain, transfer curve processing etc. + * + * A per input channel shape modification model. + * Each transfer curve adjusted value is further + * adjusted in a way that depends on the combination + * of value of all the other input channels. + * This allows for ink interaction effects in a way + * that should allow good conformance to all combinations + * of input values at 50% values. + * + * An n-linear interpolation between all combination of + * the 0 and 100% colorant primary combinations (Neugenbauer). + * + * The model making can be rather slow, particularly if high + * quality, large number of colorants, large number of sample + * points. + * + * This code is based on profile.c, sprof.c and xlut.c + * + */ + +/* + * TTBD: + * + * !!!! Should change device transfer model to include offset & scale, + * to better match display & other devices !!!! + * + * Should add Jab pcs mode, so that Jab gamuts can be + * written. + * + * Remove #ifndef DEBUG & replace with verbose progress bars. + * + * Add support for extra profile details to create() to support + * profxinf like stuff. + * + * Add ink order and overlay modeling stuff back in, with + * new ink overlay model (see mpprof0.c). + * + * Rather than computing XYZ based versios of the print model + * and ink mixing models, should compute spectrally sharpened + * equivalents to XYZ ?? (The spectral CIE base values don't + * seem much more accurate than XYZ, so perhaps remove SHARPEN code ?). + * + * Need to cleanup error handling. + */ + +#define VERSION "1.0" + +#undef DEBUG +#undef TESTDFUNC /* Check delta functions */ +#undef NODDV /* Use (very slow) non d/dv powell */ +#undef DOPLOT /* Plot the device curves */ + +#undef NOPROCESS /* Define to skip all fitting */ +#define MULTIPASS /* Fit passes by parts (normal mode) */ +#undef BIGBANG /* Define to fit all parameters at once */ +#undef ISHAPE /* define to try SVD init of shape parameters */ +#undef SHARPEN /* use sharpened XYZ values for modeling */ + /* (Wrecks derivative lookup at the moment, and makes */ + /* accuracy worse!) */ + + /* Transfer curve parameter (wiggle) minimisation weight */ +#define TRANS_BASE 0.2 /* 0 & 1 harmonic parameter weight */ +#define TRANS_HBASE 0.8 /* 2nd harmonic and above base parameter weight */ +#define SHAPE_PMW 0.2 /* Shape parameter (wiggle) minimisation weight */ +#define COMB_PMW 0.008 /* Primary combination anchor point distance weight */ + +#define verbo stdout + +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include <sys/types.h> +#include <time.h> +#if defined(__IBMC__) && defined(_M_IX86) +#include <float.h> +#endif +#include "numlib.h" +#include "cgats.h" +#include "icc.h" +#include "xicc.h" /* Spectral support */ +#include "xspect.h" /* Spectral support */ +#include "xcolorants.h" /* Known colorants support */ +#include "insttypes.h" +#include "gamut.h" +#include "mpp.h" +#ifdef DOPLOT +#include "plot.h" +#endif /* DOPLOT */ + +/* Forward declarations */ +static double bandval(mpp *p, int band, double *dev); +static double dbandval(mpp *p, double *dv, int band, double *dev); +static void forward(mpp *p, double *spec, double *Lab, double *XYZ, double *dev); +static int create(mpp *p, int verb, int quality, int display, double limit, inkmask devmask, + int spec_n, double spec_wl_short, double spec_wl_long, + double norm, instType itype, int nodp, mppcol *points); +static void compute_wb(mpp *p); +static void init_shape(mpp *p); + +/* Utilities */ + +#ifdef SHARPEN +/* Convert from XYZ to spectrally sharpened response */ +/* (Can this generate -ve values for real colors ??) */ +static void XYZ2sharp(double *a, double *b, double *c) { + double xyz[3], rgb[3]; + + xyz[0] = *a; + xyz[1] = *b; + xyz[2] = *c; + + rgb[0] = 0.8562 * xyz[0] + 0.3372 * xyz[1] - 0.1934 * xyz[2]; + rgb[1] = -0.8360 * xyz[0] + 1.8327 * xyz[1] + 0.0033 * xyz[2]; + rgb[2] = 0.0357 * xyz[0] - 0.0469 * xyz[1] + 1.0112 * xyz[2]; + + *a = rgb[0]; + *b = rgb[1]; + *c = rgb[2]; +} + +/* Convert from spectrally sharpened response to XYZ */ +static void sharp2XYZ(double *a, double *b, double *c) { + double xyz[3], rgb[3]; + + rgb[0] = *a; + rgb[1] = *b; + rgb[2] = *c; + + xyz[0] = 0.9873999149199270 * rgb[0] + - 0.1768250198556842 * rgb[1] + + 0.1894251049357572 * rgb[2]; + xyz[1] = 0.4504351090445316 * rgb[0] + + 0.4649328977527109 * rgb[1] + + 0.0846319932027575 * rgb[2]; + xyz[2] = -0.0139683251072516 * rgb[0] + + 0.0278065725014340 * rgb[1] + + 0.9861617526058175 * rgb[2]; + + *a = xyz[0]; + *b = xyz[1]; + *c = xyz[2]; +} +#endif /* SHARPEN */ + + +/* Method implimentations */ + +/* Write out the mpp to a CGATS format .mpp file */ +/* Return nz on error */ +static int write_mpp( +mpp *p, /* This */ +char *outname, /* Filename to write to */ +int dolab /* If NZ, write Lab values rather than XYZ */ +) { + int i, j, n; + time_t clk = time(0); + struct tm *tsp = localtime(&clk); + char *atm = asctime(tsp); /* Ascii time */ + cgats *ocg; /* CGATS structure */ + char *ident = icx_inkmask2char(p->imask, 1); + int nsetel = 0; + cgats_set_elem *setel; /* Array of set value elements */ + char buf[100]; + + atm[strlen(atm)-1] = '\000'; /* Remove \n from end */ + + /* Setup output cgats file */ + ocg = new_cgats(); /* Create a CGATS structure */ + ocg->add_other(ocg, "MPP"); /* our special type is Model Printer Profile */ + ocg->add_table(ocg, tt_other, 0); /* Start the first table */ + + ocg->add_kword(ocg, 0, "DESCRIPTOR", "Argyll Model Printer Profile, Colorant linearisation",NULL); + ocg->add_kword(ocg, 0, "ORIGINATOR", "Argyll mpp", NULL); + ocg->add_kword(ocg, 0, "CREATED",atm, NULL); + if (p->display) + ocg->add_kword(ocg, 0, "DEVICE_CLASS","DISPLAY", NULL); /* What sort of device this is */ + else { + ocg->add_kword(ocg, 0, "DEVICE_CLASS","OUTPUT", NULL); /* What sort of device this is */ + + /* Note what instrument the chart was read with, in case we */ + /* want to apply FWA when converting spectral data to CIE */ + ocg->add_kword(ocg, 0, "TARGET_INSTRUMENT", inst_name(p->itype) , NULL); + + sprintf(buf,"%5.1f",p->limit * 100.0); + ocg->add_kword(ocg, 0, "TOTAL_INK_LIMIT", buf, NULL); + } + + ocg->add_kword(ocg, 0, "COLOR_REP", ident, NULL); + + /* Record how many factors are used in device channel transfer curve */ + sprintf(buf,"%d",p->cord); + ocg->add_kword(ocg, 0, "TRANSFER_ORDERS", buf, NULL); + + /* Record if shaper parameters are being used */ + if (p->useshape) { + ocg->add_kword(ocg, 0, "USE_SHAPER", "YES", NULL); + } else { + ocg->add_kword(ocg, 0, "USE_SHAPER", "NO", NULL); + } + + /* Setup the table, which holds all the model parameters. */ + /* There is always a parameter per X Y Z or spectral band */ + ocg->add_field(ocg, 0, "PARAMETER", nqcs_t); + if (dolab) { + ocg->add_field(ocg, 0, "LAB_L", r_t); + ocg->add_field(ocg, 0, "LAB_A", r_t); + ocg->add_field(ocg, 0, "LAB_B", r_t); + } else { + ocg->add_field(ocg, 0, "XYZ_X", r_t); + ocg->add_field(ocg, 0, "XYZ_Y", r_t); + ocg->add_field(ocg, 0, "XYZ_Z", r_t); + } + nsetel = 1 + 3; + + /* Add fields for spectral values */ + if (p->spec_n > 0) { + + nsetel += p->spec_n; /* Spectral values */ + sprintf(buf,"%d", p->spec_n); + ocg->add_kword(ocg, 0, "SPECTRAL_BANDS",buf, NULL); + sprintf(buf,"%f", p->spec_wl_short); + ocg->add_kword(ocg, 0, "SPECTRAL_START_NM",buf, NULL); + sprintf(buf,"%f", p->spec_wl_long); + ocg->add_kword(ocg, 0, "SPECTRAL_END_NM",buf, NULL); + sprintf(buf,"%f", p->norm * 100.0); + ocg->add_kword(ocg, 0, "SPECTRAL_NORM",buf, NULL); + + /* Generate fields for spectral values */ + for (i = 0; i < p->spec_n; i++) { + int nm; + + /* Compute nearest integer wavelength */ + nm = (int)(p->spec_wl_short + ((double)i/(p->spec_n-1.0)) + * (p->spec_wl_long - p->spec_wl_short) + 0.5); + + sprintf(buf,"SPEC_%03d",nm); + ocg->add_field(ocg, 0, buf, r_t); + } + } + + if ((setel = (cgats_set_elem *)malloc(sizeof(cgats_set_elem) * nsetel)) == NULL) { + free(ident); + sprintf(p->err,"write_mpp: malloc of setel failed"); + return 1; + } + + /* Write out the transfer curve values */ + for (i = 0; i < p->n; i++) { /* each colorant */ + for (j = 0; j < p->cord; j++) { /* curve order values */ + + sprintf(buf,"t_%d_%d",i,j); + setel[0].c = buf; /* Parameter identifier */ + + for (n = 0; n < (3+p->spec_n); n++) + setel[1+n].d = p->tc[i][n][j]; + ocg->add_setarr(ocg, 0, setel); + } + } + + if (p->useshape) { + + /* Write out the shaper values */ + for (i = 0; i < p->nnn2; i++) { /* For all shaper values */ + int m = p->c2f[i].ink; + int k = p->c2f[i].comb; + + sprintf(buf,"s_%d_%d",m, k); + setel[0].c = buf; /* Parameter identifier */ + + for (n = 0; n < (3+p->spec_n); n++) + setel[1+n].d = p->shape[m][k][n]; + + ocg->add_setarr(ocg, 0, setel); + } + } + + /* Write out the colorant combination values */ + for (i = 0; i < p->nn; i++) { + + sprintf(buf,"c_%d",i); + setel[0].c = buf; /* Parameter identifier */ + + for (n = 0; n < (3+p->spec_n); n++) + setel[1+n].d = p->pc[i][n]; + +#ifdef SHARPEN + sharp2XYZ(&setel[1+0].d, &setel[1+1].d, &setel[1+2].d); +#endif + if (dolab) { + double ttt[3]; + ttt[0] = setel[1+0].d, ttt[1] = setel[1+1].d, ttt[2] = setel[1+2].d; + icmXYZ2Lab(&icmD50, ttt, ttt); + setel[1+0].d = ttt[0], setel[1+1].d = ttt[1], setel[1+2].d = ttt[2]; + } + ocg->add_setarr(ocg, 0, setel); + } + free(setel); + free(ident); + + /* Write it */ + if (ocg->write_name(ocg, outname)) { + strcpy(p->err, ocg->err); + return 1; + } + + ocg->del(ocg); /* Clean up */ + + return 0; +} + +/* Read in the mpp CGATS .mpp file */ +/* Return nz on error */ +static int read_mpp( +mpp *p, /* This */ +char *inname /* Filename to read from */ +) { + int i, j, n, ix; + cgats *icg; /* input cgats structure */ + int ti; /* Temporary CGATs index */ + int islab = 0; /* nz if Lab parameters */ + + /* Open and look at the .mpp model printer profile */ + if ((icg = new_cgats()) == NULL) { /* Create a CGATS structure */ + sprintf(p->err, "read_mpp: new_cgats() failed"); + return 2; + } + icg->add_other(icg, "MPP"); /* our special type is Model Printer Profile */ + + if (icg->read_name(icg, inname)) { + strcpy(p->err, icg->err); + icg->del(icg); + return 1; + } + + if (icg->ntables == 0 || icg->t[0].tt != tt_other || icg->t[0].oi != 0) { + sprintf(p->err, "read_mpp: Input file '%s' isn't a MPP format file",inname); + icg->del(icg); + return 1; + } + if (icg->ntables != 1) { + sprintf(p->err, "Input file '%s' doesn't contain exactly one table",inname); + icg->del(icg); + return 1; + } + if ((ti = icg->find_kword(icg, 0, "COLOR_REP")) < 0) { + sprintf(p->err, "read_mpp: Input file '%s' doesn't contain keyword COLOR_REP",inname); + icg->del(icg); + return 1; + } + + p->imask = icx_char2inkmask(icg->t[0].kdata[ti]); + p->n = icx_noofinks(p->imask); + p->nn = 1 << p->n; + p->nnn2 = p->n * p->nn/2; + + if (p->n == 0) { + sprintf(p->err, "read_mpp: COLOR_REP '%s' invalid from file '%s' (No matching devmask)", + icg->t[0].kdata[ti], inname); + icg->del(icg); + return 1; + } + + /* See if it is the expected device class */ + if ((ti = icg->find_kword(icg, 0, "DEVICE_CLASS")) < 0) { + sprintf(p->err, "read_mpp: Input file '%s' doesn't contain keyword DEVICE_CLASS",inname); + icg->del(icg); + return 1; + } + if (strcmp(icg->t[0].kdata[ti],"OUTPUT") == 0) { + + if ((ti = icg->find_kword(icg, 0, "TOTAL_INK_LIMIT")) >= 0) { + double imax; + imax = atof(icg->t[0].kdata[ti]); + p->limit = imax/100.0; + } else { + p->limit = 0.0; /* Don't use ink limit */ + } + + if ((ti = icg->find_kword(icg, 0, "TARGET_INSTRUMENT")) < 0) { + sprintf(p->err, "read_mpp: Can't find keyword TARGET_INSTRUMENT in file '%s'", inname); + icg->del(icg); + return 1; + } + + if ((p->itype = inst_enum(icg->t[0].kdata[ti])) == instUnknown + && icg->find_kword(icg, 0, "SPECTRAL_BANDS") >= 0) { + sprintf(p->err, "read_mpp: Unrecognised target instrument '%s' in file '%s'", + icg->t[0].kdata[ti], inname); + icg->del(icg); + return 1; + } + + p->display = 0; + + } else if (strcmp(icg->t[0].kdata[ti],"DISPLAY") == 0) { + + p->display = 1; + p->itype = instUnknown; + p->limit = p->n; + + } else { + /* Don't know anything else at the moment */ + sprintf(p->err, "read_mpp: Input file '%s' has unknown DEVICE_CLASS '%s'", + inname, icg->t[0].kdata[ti]); + icg->del(icg); + return 1; + } + + /* Read the number of device linearisation orders */ + if ((ti = icg->find_kword(icg, 0, "TRANSFER_ORDERS")) < 0) { + sprintf(p->err, "read_mpp: Input file '%s' doesn't contain keyword TRANSFER_ORDERS", + inname); + icg->del(icg); + return 1; + } + p->cord = atoi(icg->t[0].kdata[ti]); + if (p->cord < 1 || p->cord > MPP_MXTCORD) { + sprintf(p->err, "read_mpp: Input file '%s' has out of range TRANSFER_ORDERS %d", + inname, p->cord); + icg->del(icg); + return 1; + } + + /* See if shaper parameters are used */ + p->useshape = 0; + if ((ti = icg->find_kword(icg, 0, "USE_SHAPER")) >= 0) { + if(strcmp(icg->t[0].kdata[ti], "YES") == 0) + p->useshape = 1; + } + + /* Read the model parameters */ + { + int ci; /* Parameter dentified index */ + int spi[3+MPP_MXBANDS]; /* CGATS indexes for each band */ + char *xyzfname[3] = { "XYZ_X", "XYZ_Y", "XYZ_Z" }; + char *labfname[3] = { "LAB_L", "LAB_A", "LAB_B" }; + char buf[100]; + + /* See if we have spectral information available */ + if (icg->find_kword(icg, 0, "SPECTRAL_BANDS") < 0) { + p->spec_n = 0; /* None */ + } else { + if ((ti = icg->find_kword(icg, 0, "SPECTRAL_BANDS")) < 0) + error ("Input file doesn't contain keyword SPECTRAL_BANDS"); + p->spec_n = atoi(icg->t[0].kdata[ti]); + if ((ti = icg->find_kword(icg, 0, "SPECTRAL_START_NM")) < 0) + error ("Input file doesn't contain keyword SPECTRAL_START_NM"); + p->spec_wl_short = atof(icg->t[0].kdata[ti]); + if ((ti = icg->find_kword(icg, 0, "SPECTRAL_END_NM")) < 0) + error ("Input file doesn't contain keyword SPECTRAL_END_NM"); + p->spec_wl_long = atof(icg->t[0].kdata[ti]); + if ((ti = icg->find_kword(icg, 0, "SPECTRAL_NORM")) < 0) + error ("Input file doesn't contain keyword SPECTRAL_NORM"); + p->norm = atof(icg->t[0].kdata[ti])/100.0; + } + + if ((new_mppcol(&p->white, p->n, p->spec_n)) != 0) { + error("Malloc failed!"); + } + if ((new_mppcol(&p->black, p->n, p->spec_n)) != 0) { + error("Malloc failed!"); + } + if ((new_mppcol(&p->kblack, p->n, p->spec_n)) != 0) { + error("Malloc failed!"); + } + init_shape(p); /* Allocate and init shape related parameter space */ + + /* Get the field indexes */ + if ((ci = icg->find_field(icg, 0, "PARAMETER")) < 0) { + sprintf(p->err, "read_mpp: Input file '%s' doesn't contain field PARAMETER", + inname); + icg->del(icg); + return 1; + } + if (icg->t[0].ftype[ci] != nqcs_t) { + sprintf(p->err, "read_mpp: Input file '%s' field PARAMETER is wrong type", + inname); + icg->del(icg); + return 1; + } + + for (i = 0; i < 3; i++) { /* XYZ fields */ + if ((spi[i] = icg->find_field(icg, 0, xyzfname[i])) < 0) { + break; + } + if (icg->t[0].ftype[spi[i]] != r_t) { + sprintf(p->err, "read_mpp: Input file '%s' field %s is wrong type", + inname, buf); + icg->del(icg); + return 1; + } + } + + if (i < 3) { + islab = 1; + for (i = 0; i < 3; i++) { /* XYZ fields */ + if ((spi[i] = icg->find_field(icg, 0, labfname[i])) < 0) { + sprintf(p->err, "read_mpp: Input file '%s' doesn't contain field %s or %s", + inname, xyzfname[i], labfname[i]); + icg->del(icg); + return 1; + } + if (icg->t[0].ftype[spi[i]] != r_t) { + sprintf(p->err, "read_mpp: Input file '%s' field %s is wrong type", + inname, buf); + icg->del(icg); + return 1; + } + } + } + + /* Find the fields for spectral values */ + if (p->spec_n > 0) { + for (j = 0; j < p->spec_n; j++) { + int nm; + + /* Compute nearest integer wavelength */ + nm = (int)(p->spec_wl_short + ((double)j/(p->spec_n-1.0)) + * (p->spec_wl_long - p->spec_wl_short) + 0.5); + + sprintf(buf,"SPEC_%03d",nm); + + if ((spi[3+j] = icg->find_field(icg, 0, buf)) < 0) { + sprintf(p->err, "read_mpp: Input file '%s' doesn't contain field %s", + buf,inname); + icg->del(icg); + return 1; + } + if (icg->t[0].ftype[spi[3+j]] != r_t) { + sprintf(p->err, "read_mpp: Input file '%s' field %s is wrong type", + inname, buf); + icg->del(icg); + return 1; + } + } + } + + /* Read the transfer curve values */ + for (i = 0; i < p->n; i++) { /* each colorant */ + for (j = 0; j < p->cord; j++) { /* curve order values */ + + sprintf(buf,"t_%d_%d",i,j); + + /* Find the right parameter */ + for (ix = 0; ix < icg->t[0].nsets; ix++) { + + if (strcmp((char *)icg->t[0].fdata[ix][ci], buf) == 0) { + for (n = 0; n < (3+p->spec_n); n++) + p->tc[i][n][j] = *((double *)icg->t[0].fdata[ix][spi[n]]); + break; + } + } + } + } + + if (p->useshape) { + + /* Read the shaper values */ + for (i = 0; i < p->nnn2; i++) { /* For all shaper values */ + int m = p->c2f[i].ink; + int k = p->c2f[i].comb; + + sprintf(buf,"s_%d_%d",m, k); + + /* Find the right parameter */ + for (ix = 0; ix < icg->t[0].nsets; ix++) { + + if (strcmp((char *)icg->t[0].fdata[ix][ci], buf) == 0) { + for (n = 0; n < (3+p->spec_n); n++) + p->shape[m][k][n] = *((double *)icg->t[0].fdata[ix][spi[n]]); + + break; + } + } + } + } + + /* Read the combination values */ + for (i = 0; i < p->nn; i++) { + + sprintf(buf,"c_%d",i); + + + /* Find the right parameter */ + for (ix = 0; ix < icg->t[0].nsets; ix++) { + + if (strcmp((char *)icg->t[0].fdata[ix][ci], buf) == 0) { + for (n = 0; n < (3+p->spec_n); n++) + p->pc[i][n] = *((double *)icg->t[0].fdata[ix][spi[n]]); + if (islab) { + double tt[3]; + tt[0] = p->pc[i][0], tt[1] = p->pc[i][1], tt[2] = p->pc[i][2]; + icmLab2XYZ(&icmD50, tt, tt); + p->pc[i][0] = tt[0], p->pc[i][1] = tt[1], p->pc[i][2] = tt[2]; + } +#ifdef SHARPEN + XYZ2sharp(&p->pc[i][0], &p->pc[i][1], &p->pc[i][2]); +#endif + break; + } + } + } + } + + icg->del(icg); /* Clean up */ + + /* Figure out the white and black points */ + compute_wb(p); + + return 0; +} + +/* Get various types of information about the mpp */ +static void get_info( +mpp *p, /* This */ +inkmask *imask, /* Inkmask, describing device colorspace */ +int *nodchan, /* Number of device channels */ +double *limit, /* Total ink limit (0.0 .. devchan) */ +int *spec_n, /* Number of spectral bands, 0 if none */ +double *spec_wl_short, /* First reading wavelength in nm (shortest) */ +double *spec_wl_long, /* Last reading wavelength in nm (longest) */ +instType *itype, /* Instrument type */ +int *display /* Return nz if display type */ +) { + if (imask != NULL) + *imask = p->imask; + if (nodchan != NULL) + *nodchan = p->n; + if (limit != NULL) + *limit = p->limit; + if (spec_n != NULL) + *spec_n = p->spec_n; + if (spec_wl_short != NULL) + *spec_wl_short = p->spec_wl_short; + if (spec_wl_long != NULL) + *spec_wl_long = p->spec_wl_long; + if (itype != NULL) + *itype = p->itype; + if (display != NULL) + *display = p->display; +} + +/* Set an illuminant and observer to use spectral model */ +/* for CIE lookup with optional FWA. Set both to default for XYZ mpp model. */ +/* return 0 on OK, 1 on spectral not supported, 2 on other error */ +/* If the model is for a display, the illuminant will be ignored. */ +static int set_ilob( +mpp *p, +icxIllumeType ilType, /* Illuminant type (icxIT_default for none) */ +xspect *custIllum, /* Custom illuminant (NULL for none) */ +icxObserverType obType, /* Observer type (icxOT_default for none) */ +xspect custObserver[3], /* Custom observer (NULL for none) */ +icColorSpaceSignature rcs, /* Return color space, icSigXYZData or icSigLabData */ +int use_fwa /* NZ to involke FWA. */ +) { + + /* Get rid of any existing conversion object */ + if (p->spc != NULL) { + p->spc->del(p->spc); + p->spc = NULL; + } + + p->pcs = rcs; + + if (ilType == icxIT_default && obType == icxOT_default && use_fwa == 0) + return 0; + + if (p->spec_n == 0) { + p->errc = 1; + sprintf(p->err,"No Spectral Data in MPP"); + return 1; + } + + if (p->display) { + ilType = icxIT_none; + custIllum = NULL; + } + + if ((p->spc = new_xsp2cie(ilType, custIllum, obType, custObserver, rcs, 1)) == NULL) + error("mpp->set_ilob, new_xsp2cie failed"); + + if (use_fwa) { + int j; + xspect white, inst; + + white.norm = p->norm; + white.spec_n = p->spec_n; + white.spec_wl_short = p->spec_wl_short; + white.spec_wl_long = p->spec_wl_long; + for (j = 0; j < p->spec_n; j++) + white.spec[j] = p->white.band[3+j]; + + if (inst_illuminant(&inst, p->itype) != 0) + error ("mpp->set_ilob, instrument doesn't have an FWA illuminent"); + + if (p->spc->set_fwa(p->spc, &inst, NULL, &white)) + error ("mpp->set_ilob, set_fwa faild"); + } + + return 0; +} + +/* Lookup an XYZ or Lab color */ +static void lookup( +mpp *p, /* This */ +double *out, /* Returned XYZ or Lab */ +double *in /* Input device values */ +) { + if (p->spc == NULL) { + double *Lab; + double *XYZ; + + if (p->pcs == icSigLabData) { + Lab = out; + XYZ = NULL; + } else { + Lab = NULL; + XYZ = out; + } + + forward(p, NULL, Lab, XYZ, in); + + return; + } else { + xspect tspec; + + p->lookup_spec(p, &tspec, in); + p->spc->convert(p->spc, out, &tspec); + } +} + +/* Lookup an XYZ or Lab color, plus the partial derivative. */ +/* This is useful if the lookup is being used within */ +/* a minimisation routine */ +static void dlookup( +mpp *p, +double *out, /* Return the XYZ or Lab */ +double **dout, /* Return the partial derivative dout[3][n] */ +double *dev) { + int i, j, k; + +#ifdef SHARPEN + /* Can't handle this without converting from sharpened derivative */ + /* to XYZ derivative. ~~~~~~ */ +#endif + /* We can't handle derivative using FWA, so */ + /* always compute from the XYZ model values. */ + for (j = 0; j < 3; j++) { /* Compute each bands value */ + out[j] = dbandval(p, dout[j], j, dev); + } + if (p->pcs == icSigLabData) { + double dlab[3][3]; + + icxdXYZ2Lab(&icmD50, out, dlab, out); + + /* Apply Lab deriv to dout */ + for (i = 0; i < p->n; i++) { + double tt[3]; + + for (k = 0; k < 3; k++) + tt[k] = dout[k][i]; + + for (j = 0; j < 3; j++) { + dout[j][i] = 0.0; + + for (k = 0; k < 3; k++) { + dout[j][i] += dlab[j][k] * tt[k]; + } + } + } + } +} + + +/* Get the white and black points */ +static void get_wb( +mpp *p, /* This */ +double *white, +double *black, +double *kblack /* K only black */ +) { + if (white != NULL) { + if (p->spc == NULL) { + white[0] = p->white.band[0]; + white[1] = p->white.band[1]; + white[2] = p->white.band[2]; +#ifdef SHARPEN + sharp2XYZ(&white[0], &white[1], &white[2]); +#endif + if (p->pcs == icSigLabData) + icmXYZ2Lab(&icmD50, white, white); + } else { + int j; + xspect ispect; + + ispect.norm = p->norm; + ispect.spec_n = p->spec_n; + ispect.spec_wl_short = p->spec_wl_short; + ispect.spec_wl_long = p->spec_wl_long; + for (j = 0; j < p->spec_n; j++) + ispect.spec[j] = p->white.band[3+j]; + + p->spc->convert(p->spc, white, &ispect); + } + } + if (black != NULL) { + if (p->spc == NULL) { + black[0] = p->black.band[0]; + black[1] = p->black.band[1]; + black[2] = p->black.band[2]; +#ifdef SHARPEN + sharp2XYZ(&black[0], &black[1], &black[2]); +#endif + if (p->pcs == icSigLabData) + icmXYZ2Lab(&icmD50, black, black); + } else { + int j; + xspect ispect; + + ispect.norm = p->norm; + ispect.spec_n = p->spec_n; + ispect.spec_wl_short = p->spec_wl_short; + ispect.spec_wl_long = p->spec_wl_long; + for (j = 0; j < p->spec_n; j++) + ispect.spec[j] = p->black.band[3+j]; + + p->spc->convert(p->spc, black, &ispect); + } + } + + if (kblack != NULL) { + if (p->spc == NULL) { + kblack[0] = p->kblack.band[0]; + kblack[1] = p->kblack.band[1]; + kblack[2] = p->kblack.band[2]; +#ifdef SHARPEN + sharp2XYZ(&kblack[0], &kblack[1], &kblack[2]); +#endif + if (p->pcs == icSigLabData) + icmXYZ2Lab(&icmD50, kblack, kblack); + } else { + int j; + xspect ispect; + + ispect.norm = p->norm; + ispect.spec_n = p->spec_n; + ispect.spec_wl_short = p->spec_wl_short; + ispect.spec_wl_long = p->spec_wl_long; + for (j = 0; j < p->spec_n; j++) + ispect.spec[j] = p->kblack.band[3+j]; + + p->spc->convert(p->spc, kblack, &ispect); + } + } +} + +/* Lookup an XYZ value. */ +/* (Note that this is never FWA compensated) */ +static void lookup_xyz( +mpp *p, /* This */ +double *out, /* Returned XYZ value */ +double *in /* Input device values */ +) { + forward(p, NULL, NULL, out, in); +} + +/* Lookup a spectral value. */ +/* (Note that this is never FWA compensated) */ +static void lookup_spec( +mpp *p, /* This */ +xspect *out, /* Returned spectral value */ +double *in /* Input device values */ +) { + int j; + + out->norm = p->norm; + out->spec_n = p->spec_n; + out->spec_wl_short = p->spec_wl_short; + out->spec_wl_long = p->spec_wl_long; + + forward(p, out->spec, NULL, NULL, in); + + for (j = 0; j < p->spec_n; j++) + out->spec[j] *= out->norm; +} + +/* Macros for an arbitrary dimensional counter */ +/* Declare the counter name nn, dimensions di, & count */ + +#define DCOUNT(nn, di, start, reset, count) \ + int nn[MAX_CHAN]; /* counter value */ \ + int nn##_di = (di); /* Number of dimensions */ \ + int nn##_stt = (start); /* start count value */ \ + int nn##_rst = (reset); /* reset on carry value */ \ + int nn##_res = (count); /* last count +1 */ \ + int nn##_e /* dimension index */ + +/* Set the counter value to 0 */ +#define DC_INIT(nn) \ +{ \ + for (nn##_e = 0; nn##_e < nn##_di; nn##_e++) \ + nn[nn##_e] = nn##_stt; \ + nn##_e = 0; \ +} + +/* Increment the counter value */ +#define DC_INC(nn) \ +{ \ + for (nn##_e = 0; nn##_e < nn##_di; nn##_e++) { \ + nn[nn##_e]++; \ + if (nn[nn##_e] < nn##_res) \ + break; /* No carry */ \ + nn[nn##_e] = nn##_rst; \ + } \ +} + +/* After increment, expression is TRUE if counter is done */ +#define DC_DONE(nn) \ + (nn##_e >= nn##_di) + +/* Create a gamut */ +/* (This will be in the current PCS) */ +static gamut *mpp_gamut( +mpp *p, /* This */ +double detail /* gamut detail level, 0.0 = def */ +) { + gamut *gam; + int res; /* Sample point resolution */ + double white[3], black[3], kblack[3]; + DCOUNT(co, p->n, 0, 0, 2); + + if (detail == 0.0) + detail = 10.0; + + gam = new_gamut(detail, 0, 0); /* Lab only at the moment */ + + /* Explore the gamut by itterating through */ + /* it with sample points in device space. */ + + res = (int)(100.0/detail); /* Establish an appropriate sampling density */ + + if (res < 3) + res = 3; + + DC_INIT(co); + + /* Itterate over all the 2 faces in the device space. */ + /* We're assuming in practice that only colors lying */ + /* on such faces contribute to the gamut volume. */ + /* The total number of faces explored will be */ + /* (dim * (dim-1))/2 * 2 ^ (dim-2) */ + /* It seems possible that a number of these faces could */ + /* be cheaply culled by examining the PCS corner values */ + /* for each 2^(dim-2) combinations, and skiping those. */ + /* that cannot possibly expand the gamut ? */ + /* ie. that they are clearly contained by other face combinations. */ + while(!DC_DONE(co)) { + int e, m1, m2; + double in[MAX_CHAN]; + double out[3]; + double sum; + + /* Check the ink limit */ + for (sum = 0.0, e = 0; e < p->n; e++) + sum += (double)co[e]; + + if (p->limit > 1e-4 && (sum - 1.0) > p->limit) { + DC_INC(co); + continue; /* Skip points really over limit */ + } + + /* Scan only device surface */ + for (m1 = 0; m1 < p->n; m1++) { + if (co[m1] != 0) + continue; + for (m2 = m1 + 1; m2 < p->n; m2++) { + int x, y; + + if (co[m2] != 0) + continue; + + for (e = 0; e < p->n; e++) + in[e] = (double)co[e]; /* Base value */ + + for (x = 0; x < res; x++) { /* step over surface */ + in[m1] = x/(res - 1.0); + for (y = 0; y < res; y++) { + double ssum, iin[MAX_CHAN]; + in[m2] = y/(res - 1.0); + ssum = sum + in[m1] + in[m2]; + if (p->limit > 1e-4 && (ssum - 1.0) > p->limit) { + continue; + } + + for (e = 0; e < p->n; e++) + iin[e] = in[e]; /* Scalable copy */ + + /* Apply ink limit by simple scaling */ + if (p->limit > 1e-4 && ssum > p->limit) { + for (e = 0; e < p->n; e++) + iin[e] *= p->limit/ssum; + } + + forward(p, NULL, out, NULL, iin); + gam->expand(gam, out); + } + } + } + } + + /* Increment index within block */ + DC_INC(co); + } + + /* set the white and black points */ + p->get_wb(p, white, black, kblack); + gam->setwb(gam, white, black, kblack); + + /* Set the cusp points */ + /* Do this by scanning just 0 & 100% colorant combinations */ + + res = 2; + gam->setcusps(gam, 0, NULL); + DC_INIT(co); + + /* Itterate over all the faces in the device space */ + while(!DC_DONE(co)) { + int e, m1, m2; + double in[MAX_CHAN]; + double out[3]; + double sum; + + /* Check the ink limit */ + for (sum = 0.0, e = 0; e < p->n; e++) + sum += (double)co[e]; + + if (p->limit > 1e-4 && (sum - 1.0) > p->limit) { + DC_INC(co); + continue; /* Skip points really over limit */ + } + + /* Scan only device surface */ + for (m1 = 0; m1 < p->n; m1++) { + if (co[m1] != 0) + continue; + for (m2 = m1 + 1; m2 < p->n; m2++) { + int x, y; + + if (co[m2] != 0) + continue; + + for (e = 0; e < p->n; e++) + in[e] = (double)co[e]; /* Base value */ + + for (x = 0; x < res; x++) { /* step over surface */ + in[m1] = x/(res - 1.0); + for (y = 0; y < res; y++) { + double ssum, iin[MAX_CHAN]; + in[m2] = y/(res - 1.0); + ssum = sum + in[m1] + in[m2]; + if (p->limit > 1e-4 && (ssum - 1.0) > p->limit) { + continue; + } + + for (e = 0; e < p->n; e++) + iin[e] = in[e]; /* Scalable copy */ + + /* Apply ink limit by simple scaling */ + if (p->limit > 1e-4 && ssum > p->limit) { + for (e = 0; e < p->n; e++) + iin[e] *= p->limit/ssum; + } + + forward(p, NULL, out, NULL, iin); + gam->setcusps(gam, 1, out); + } + } + } + } + + /* Increment index within block */ + DC_INC(co); + } + gam->setcusps(gam, 2, NULL); + + return gam; +} + +/* Delete it */ +static void del_mpp(mpp *p) { + if (p != NULL) { + del_mppcol(&p->white, p->n, p->spec_n); + del_mppcol(&p->black, p->n, p->spec_n); + del_mppcol(&p->kblack, p->n, p->spec_n); + del_mppcols(p->cols, p->nodp, p->n, p->spec_n); /* Delete array of target points */ + if (p->spc != NULL) + p->spc->del(p->spc); + + /* Delete shape parameters */ + if (p->shape != NULL) { + int i, j; + for (j = 0; j < p->n; j++) { + if (p->shape[j] != NULL) { + for (i = 0; i < p->nn; i++) { + if (p->shape[j][i] != NULL) + free(p->shape[j][i]); + } + free(p->shape[j]); + } + } + free(p->shape); + } + free(p); + } +} + +/* Allocate a new, uninitialised mpp */ +/* Note thate black and white points aren't allocated */ +mpp *new_mpp(void) { + mpp *p; + + if ((p = (mpp *)calloc(1, sizeof(mpp))) == NULL) + return NULL; + + p->pcs = icSigXYZData; + + /* Init method pointers */ + p->del = del_mpp; + p->create = create; + p->get_gamut = mpp_gamut; + p->write_mpp = write_mpp; + p->read_mpp = read_mpp; + p->get_info = get_info; + p->set_ilob = set_ilob; + p->get_wb = get_wb; + p->lookup = lookup; + p->dlookup = dlookup; + p->lookup_xyz = lookup_xyz; + p->lookup_spec = lookup_spec; + + return p; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Utility functions */ + +/* Allocate the data portion of an mppcol */ +/* Return NZ if malloc error */ +int new_mppcol( +mppcol *p, /* mppcol to allocate */ +int n, /* Number of inks */ +int nb /* Number of spectral bands */ +) { + int nn = (1 << n); /* Number of ink combinations */ + + if ((p->nv = (double *)malloc(sizeof(double) * n)) == NULL) { + del_mppcol(p, n, nb); + return 1; + } + if ((p->band = (double *)malloc(sizeof(double) * (3 + nb))) == NULL) { + del_mppcol(p, n, nb); + return 1; + } + if ((p->lband = (double *)malloc(sizeof(double) * (3 + nb))) == NULL) { + del_mppcol(p, n, nb); + return 1; + } + if ((p->tcnv = (double *)calloc(n, sizeof(double))) == NULL) { + del_mppcol(p, n, nb); + return 1; + } + if ((p->scnv = (double *)calloc(n, sizeof(double))) == NULL) { + del_mppcol(p, n, nb); + return 1; + } + if ((p->pcnv = (double *)malloc(nn * sizeof(double))) == NULL) { + del_mppcol(p, n, nb); + return 1; + } + if ((p->fcnv = (double *)malloc(n * nn/2 * sizeof(double))) == NULL) { + del_mppcol(p, n, nb); + return 1; + } + return 0; +} + +/* Copy the contents of one mppcol to another */ +void copy_mppcol( +mppcol *d, /* Destination */ +mppcol *s, /* Source */ +int n, /* Number of inks */ +int nb /* Number of spectral bands */ +) { + mppcol al; + int nn = (1 << n); /* Number of ink combinations */ + int i, j; + + al = *d; /* Copy destinations allocations */ + *d = *s; /* Copy source contents & pointers */ + + /* Now fixup allocation addresses */ + d->nv = al.nv; + d->band = al.band; + d->lband = al.lband; + d->tcnv = al.tcnv; + d->scnv = al.scnv; + d->pcnv = al.pcnv; + d->fcnv = al.fcnv; + + /* Copy contents of allocated data */ + for (j = 0; j < n; j++) + d->nv[j] = s->nv[j]; + + for (j = 0; j < (3+nb); j++) + d->band[j] = s->band[j]; + + for (j = 0; j < (3+nb); j++) + d->lband[j] = s->lband[j]; + + for (i = 0; i < n; i++) + d->tcnv[i] = s->tcnv[i]; + + for (i = 0; i < n; i++) + d->scnv[i] = s->scnv[i]; + + for (i = 0; i < nn; i++) + d->pcnv[i] = s->pcnv[i]; + + for (i = 0; i < (n * nn/2); i++) + d->fcnv[i] = s->fcnv[i]; +} + + +/* Free the data allocation of an mppcol */ +void del_mppcol( +mppcol *p, +int n, /* Number of inks */ +int nb /* Number of spectral bands */ +) { + if (p != NULL) { + if (p->nv != NULL) + free (p->nv); + + if (p->band != NULL) + free (p->band); + + if (p->lband != NULL) + free (p->lband); + + if (p->tcnv != NULL) + free (p->tcnv); + + if (p->scnv != NULL) + free (p->scnv); + + if (p->pcnv != NULL) + free (p->pcnv); + + if (p->fcnv != NULL) + free (p->fcnv); + } +} + +/* Allocate an array of mppcol */ +/* Return NULL if malloc error */ +mppcol *new_mppcols( +int no, /* Number in array */ +int n, /* Number of inks */ +int nb /* Number of spectral bands */ +) { + mppcol *p; + int i; + + if ((p = (mppcol *)calloc(no, sizeof(mppcol))) == NULL) { + return NULL; + } + + for (i = 0; i < no; i++) { + if (new_mppcol(&p[i], n, nb) != 0) { + del_mppcols(p, no, n, nb); + return NULL; + } + } + + return p; +} + +/* Free an array of mppcol */ +void del_mppcols( +mppcol *p, +int no, /* Number in array */ +int n, /* Number of inks */ +int nb /* Number of spectral bands */ +) { + if (p != NULL) { + int i; + + for (i = 0; i < no; i++) + del_mppcol(&p[i], n, nb); + free (p); + } +} + +/* Allocate and setup the extra shape parameters combinations */ +static void init_shape(mpp *p) { + int i, j, k; + int ix[MPP_MXINKS]; + + /* First allocate the shape parameter array */ + + if ((p->shape = (double ***)malloc(p->n * sizeof(double **))) == NULL) + error("Malloc failed (mpp shape)!"); + for (j = 0; j < p->n; j++) { + if ((p->shape[j] = (double **)malloc(p->nn * sizeof(double *))) == NULL) + error("Malloc failed (mpp shape)!"); + for (i = 0; i < p->nn; i++) { + if ((i & (1 << j)) == 0) { /* Valid combo for this ink */ + if ((p->shape[j][i] = (double *)malloc((3+p->spec_n) * sizeof(double))) == NULL) + error("Malloc failed (mpp shape)!"); + for (k = 0; k < (3+p->spec_n); k++) + p->shape[j][i][k] = 0.0; /* Initial shape value */ + } else { + p->shape[j][i] = NULL; + } + } + } + + /* Setup sparse to full and back indexing lookup */ + for (j = 0; j < p->n; j++) + ix[j] = 0; + + for (i = 0; i < p->nn; i++) { + for (j = 0; j < p->n; j++) { + p->f2c[j][i] = j * p->nn/2 + ix[j]; + if ((i & (1 << j)) == 0) { + p->c2f[j * p->nn/2 + ix[j]].ink = j; + p->c2f[j * p->nn/2 + ix[j]].comb = i; + ix[j]++; + } + } + } + +#ifdef NEVER + /* Print result */ + for (j = 0; j < p->n; j++) { + for (i = 0; i < p->nn; i++) { + printf("f2c[%d][%d] = %d\n", j,i,p->f2c[j][i]); + } + for (i = 0; i < p->nn/2; i++) { + printf("c2f[%d].ink,comb = %d, %d\n", + j * p->nn/2 + i,p->c2f[j * p->nn/2 + i].ink, + p->c2f[j * p->nn/2 + i].comb); + } + } +#endif /* NEVER */ +} + + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Creation helper functions */ + +/* Per band error metric */ + +/* Convert argument to L* like scale */ +#define lDE(aa) ((aa) > 0.008856451586 ? 116.0 * pow((aa),1.0/3.0) - 16.0 : (aa) * 903.2962896) + +/* Function to approximate visible difference in an XYZ or spectral difference */ +static double DEsq(double aa, double bb) { + double dd; + + /* Convert input to L* like value */ + aa = lDE(aa); + bb = lDE(bb); + + dd = aa - bb; + + return dd * dd; +} + +/* Convert argument to the derivative of the DEsq function at that value */ +#define dDE(aa) ((aa) > 0.008856451586 ? 38.666667 * pow((aa), -2.0/3.0) : 903.2962896) + + +/* Given a device value, return the given bands value */ +/* according to the current model. */ +static double bandval(mpp *p, int band, double *dev) { + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int m, k; + double ov; + int j = band; + + /* Compute the tranfer corrected device values */ + for (m = 0; m < p->n; m++) { + tcnv[m] = icxTransFunc(p->tc[m][j], p->cord, dev[m]); + tcnv1[m] = 1.0 - tcnv[m]; + } + + if (p->useshape) { + + for (m = 0; m < p->n; m++) + ww[m] = 0.0; + + /* Lookup the shape values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + for (m = 0; m < p->n; m++) { + ww[m] += p->shape[m][k & ~(1<<m)][j] * vv; + /* Apply weighting to shape vertex value */ + } + } + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double vv = tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + tcnv[m] = vv; + tcnv1[m] = 1.0 - vv; + } + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double vv = p->pc[k][j]; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + return ov; +} + +/* Given a device value, return the given bands value */ +/* according to the current model, as well as the partial */ +/* derivative wrt the device values. */ +static double dbandval(mpp *p, double *dv, int band, double *dev) { + + double dtcnv_ddev[MPP_MXINKS]; /* Del in tcnv[m] due to del dev[m] */ + double dww_dtcnv[MPP_MXINKS][MPP_MXINKS]; /* Del in ww[m] due to del in tcnv[m] */ + double dtcnv_tc[MPP_MXINKS]; /* Del in tcnv'[m] due to del in tcnv[m] */ + double dtcnv_ww[MPP_MXINKS]; /* Del in tcnv'[m] due to del in ww[m] */ + double dov[MPP_MXINKS]; /* Del of ov due to del in tcnv'[m] */ + + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int m, k; + double ov; + int j = band; + + /* Compute the tranfer corrected device values */ + for (m = 0; m < p->n; m++) { + tcnv[m] = icxdiTransFunc(p->tc[m][j], &dtcnv_ddev[m], p->cord, dev[m]); + tcnv1[m] = 1.0 - tcnv[m]; + ww[m] = 0.0; + } + + /* Lookup the shape values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + for (m = 0; m < p->n; m++) { + ww[m] += p->shape[m][k & ~(1<<m)][j] * vv; + /* Apply weighting to shape vertex value */ + } + } + + /* Compute del ww[m][m] for del tcnv[m] */ + for (m = 0; m < p->n; m++) { /* For each input channel were doing del of */ + int mo; + for (mo = 0; mo < p->n; mo++) + dww_dtcnv[mo][m] = 0.0; + for (k = 0; k < p->nn; k++) { + int mm; + double vv = 1.0; + for (mm = 0; mm < p->n; mm++) { /* Compute weight for node k */ + if (m == mm) + continue; + if (k & (1 << mm)) + vv *= tcnv[mm]; + else + vv *= tcnv1[mm]; + } + for (mo = 0; mo < p->n; mo++) { /* For each output channel */ + double vvv = vv * p->shape[mo][k & ~(1<<mo)][j]; + if (k & (1 << m)) + dww_dtcnv[mo][m] += vvv; + else + dww_dtcnv[mo][m] -= vvv; + } + } + } + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double tt; + double gg = ww[m]; /* Curve adjustment */ + double sv, dsv, vv = tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + tt = gg - gg * vv + 1.0; + sv = vv/tt; + dsv = (gg + 1.0)/(tt * tt); + } else { + tt = 1.0 - gg * vv; + sv = (vv - gg * vv)/tt; + dsv = (1.0 - gg)/(tt * tt); + } + tcnv[m] = sv; + tcnv1[m] = 1.0 - sv; + dtcnv_tc[m] = dsv; /* del in tcnv[m] due to del in tcnv[m] */ + dtcnv_ww[m] = (vv * vv - vv)/(tt * tt); /* del in tcnv[m] due to del in ww[m] */ + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double vv = p->pc[k][j]; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + /* Compute del ov for del tcnv[m] */ + for (m = 0; m < p->n; m++) { + for (dov[m] = 0.0, k = 0; k < p->nn; k++) { + int mm; + double vv = p->pc[k][j]; + for (mm = 0; mm < p->n; mm++) { + if (m == mm) + continue; + if (k & (1 << mm)) + vv *= tcnv[mm]; + else + vv *= tcnv1[mm]; + } + if (k & (1 << m)) + dov[m] += vv; + else + dov[m] -= vv; + } + } + + /* Accumulate delta from input value to output */ + for (m = 0; m < p->n; m++) { + double ttt; + int mm; + + /* delta via dww */ + for (ttt = 0.0, mm = 0; mm < p->n; mm++) + ttt += dov[mm] * dtcnv_ww[mm] * dww_dtcnv[mm][m] * dtcnv_ddev[m]; + + /* delta direct */ + dv[m] = ttt + dov[m] * dtcnv_tc[m] * dtcnv_ddev[m]; + } + + return ov; +} + +/* Given a device value, return the XYZ, Lab and spectral */ +/* values according to the current model. */ +/* Pointer may be NULL if result is not needed */ +static void forward(mpp *p, double *spec, double *Lab, double *XYZ, double *dev) { + double tXYZ[3]; + int sb = 3, eb = 3; /* Start and end bands to compute */ + int j; + + if (XYZ != NULL || Lab != NULL) + sb = 0; + if (spec != NULL) + eb = 3 + p->spec_n; + + for (j = sb; j < eb; j++) { /* Compute each bands value */ + double ov; + + ov = bandval(p, j, dev); + + if (j < 3) + tXYZ[j] = ov; + else + spec[j-3] = ov; + } + +#ifdef SHARPEN + if (sb == 0) + sharp2XYZ(&tXYZ[0], &tXYZ[1], &tXYZ[2]); +#endif + if (XYZ != NULL) { + XYZ[0] = tXYZ[0]; + XYZ[1] = tXYZ[1]; + XYZ[2] = tXYZ[2]; + } + if (Lab != NULL) { + icmXYZ2Lab(&icmD50, Lab, tXYZ); /* Convert D50 Lab */ + } +} + +/* Return the specified bands current average error */ +static void banderr( +mpp *p, +double *avese, /* Return average Error from spectral bands */ +double *maxse, /* Return maximum Error from spectral bands */ +int band +) { + double serr = 0.0, smax = 0.0; + int i; + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + double spv; + mppcol *c = &p->cols[i]; + + /* Compute models output */ + spv = bandval(p, band, c->nv); + spv = sqrt(DEsq(c->band[band], spv)); + + if (spv > smax) + smax = spv; + serr += spv; + } + + if (avese != NULL) + *avese = serr/((double)p->nodp); + if (maxse != NULL) + *maxse = smax; +} + +/* Figure out what the current model errors are and return them */ +static void deltae( +mpp *p, +double *avede, /* Return average Delta E from XYZ bands */ +double *maxde, /* Return maximum Delta E from XYZ bands */ +double *avese, /* Return average Error from spectral bands */ +double *maxse /* Return maximum Error from spectral bands */ +) { + double spec[MPP_MXBANDS]; /* spectral value for each point */ + double Lab[3]; + double err = 0.0, max = 0.0; /* Delta E */ + double serr = 0.0, smax = 0.0; /* Delta S */ + int i, j; + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + double ee; + mppcol *c = &p->cols[i]; + + /* Compute models output */ + forward(p, p->spec_n > 0 ? spec : NULL, Lab, c->cXYZ, c->nv); + + /* Track the Delta E */ + c->err = icmLabDEsq(Lab, c->Lab); + ee = sqrt(c->err); + if (ee > max) + max = ee; + err += ee; + + /* Track the spectral error */ + for (j = 0; j < p->spec_n; j++) { + ee = DEsq(c->band[3+j], spec[j]); + ee = sqrt(ee); + if (ee > smax) + smax = ee; + serr += ee; + } + } + + if (avede != NULL) + *avede = err/((double)p->nodp); + if (maxde != NULL) + *maxde = max; + + if (p->spec_n > 0 && avese != NULL) + *avese = serr/(((double)p->nodp) * ((double)p->spec_n)); + if (maxse != NULL) + *maxse = smax; +} + +/* - - - - - - - - - - - - - - - */ +/* Powell optimisation callbacks */ + +/* Progress reporter for mpp powell/conjgrad */ +static void mppprog(void *pdata, int perc) { + mpp *p = (mpp *)pdata; + + if (p->verb) { + printf("%c% 3d%%",cr_char,perc); + if (perc == 100) + printf("\n"); + fflush(stdout); + } +} + +#ifdef NEVER // Skip efunc1 passes for now. + +/* Setup test point data ready for efunc1 on the given och and oba */ +static void sfunc1(mpp *p) { + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + int k = p->och; /* Channel being optimised */ + int j = p->oba; /* Band being optimised */ + int i, m; + + /* Pre-compute combination weighting and colorant combination values */ + /* for this band, without this device channels contribution. */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + int kk; + + /* Compute the tranfer corrected device values */ + for (m = 0; m < p->n; m++) { + tcnv[m] = icxTransFunc(p->tc[m][j], p->cord, c->nv[m]); + tcnv1[m] = 1.0 - tcnv[m]; + } + + /* Compute the test point primary combination values */ + c->tpcnv = c->tpcnv1 = 0.0; + for (kk = 0; kk < p->nn; kk++) { /* Combination value */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { + if (m == k) + continue; /* Multiply this in efunc1 */ + if (kk & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + if (kk & (1 << k)) + c->tpcnv += vv * p->pc[kk][j]; + else + c->tpcnv1 += vv * p->pc[kk][j]; + } + } +} + +#endif /* NEVER */ + +#ifdef NEVER // Skip efunc1 passes for now. +/* Optimise a device transfer curve to minimise a particular bands error */ +/* (Assume no shape parameters) */ +static double efunc1(void *adata, double pv[]) { + mpp *p = (mpp *)adata; + double smv, rv = 0.0; + double tcnv; /* Transfer curve corrected device values */ + int k = p->och; /* Channel being optimised */ + int j = p->oba; /* Band being optimised */ + int pc, m, i; + + /* For each test point */ + for (pc = m = 0; m < p->nodp; m++) { + mppcol *c = &p->cols[m]; + double vv; + + if (c->w < 1e-6) + continue; /* Skip zero weighted points */ + + /* Compute the tranfer corrected device values */ + tcnv = icxTransFunc(pv, p->cord, c->nv[k]); + + /* Compute band value */ + vv = tcnv * c->tpcnv + (1.0 - tcnv) * c->tpcnv1; + + /* Compute the band value error */ + vv = lDE(vv) - c->lband[j]; + rv += c->w * vv * vv; + pc++; + } + rv /= (double)pc; + + /* Compute average magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + for (smv = 0.0, m = 0; m < p->cord; m++) { + smv += pv[m] * pv[m]; + } + smv /= (double)(p->cord); + rv += 0.01 * smv; + +#ifdef DEBUG + printf("efunc1 itt %d/%d chan %d band %d k0 %f returning %f\n",p->oit,p->ott,k,j,pv[0],rv); +#endif + return rv; +} +#endif /* NEVER */ + +/* Optimise all transfer curves simultaniously to minimise a particular bands error */ +static double efunc2(void *adata, double pv[]) { + mpp *p = (mpp *)adata; + double smv, rv = 0.0; + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int j = p->oba; /* Band being optimised */ + int i, m, k; + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double ov; + + /* Compute the tranfer corrected device values */ + for (m = 0; m < p->n; m++) { + tcnv[m] = icxTransFunc(&pv[m * p->cord], p->cord, c->nv[m]); + tcnv1[m] = 1.0 - tcnv[m]; + } + + if (p->useshape) { + + for (m = 0; m < p->n; m++) + ww[m] = 0.0; + + /* Lookup the shape values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + for (m = 0; m < p->n; m++) { + ww[m] += p->shape[m][k & ~(1<<m)][j] * vv; + /* Apply weighting to shape vertex value */ + } + } + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double vv = tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + tcnv[m] = vv; + tcnv1[m] = 1.0 - vv; + } + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double vv = p->pc[k][j]; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + /* Compute the band value error */ + ov = lDE(ov) - c->lband[j]; + rv += ov * ov; + } + + rv /= (double)p->nodp; + + /* Compute weighted magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + for (smv = 0.0, m = 0; m < p->n; m++) { + double w; + for (k = 0; k < p->cord; k++) { + i = m * p->cord + k; + w = (k < 2) ? TRANS_BASE : k * TRANS_HBASE; /* Increase weight with harmonics */ + smv += w * pv[i] * pv[i]; + } + } + smv /= (double)(p->n); + rv += smv; + +#ifdef DEBUG + printf("efunc2 itt %d/%d band %d returning %f\n",p->oit,p->ott,j,rv); +#endif + return rv; +} + +/* Return the gradient of the minimisation function at the given location, */ +/* as well as the function value at this location. */ +static double dfunc2(void *adata, double dv[], double pv[]) { + mpp *p = (mpp *)adata; + double smv, tt, rv = 0.0; + double dtcnv_dpv[MPP_MXINKS][MPP_MXTCORD]; /* Del in tcnv[m] due to del in parameter */ + double dww_dtcnv[MPP_MXINKS][MPP_MXINKS]; /* Del in ww[m] due to del in tcnv[m] */ + double dtcnv_tc[MPP_MXINKS]; /* Del in tcnv'[m] due to del in tcnv[m] */ + double dtcnv_ww[MPP_MXINKS]; /* Del in tcnv'[m] due to del in ww[m] */ + double dov[MPP_MXINKS]; /* Del of ov due to del in tcnv'[m] */ + double ddov; /* Del in final ov due to del in raw ov */ + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int j = p->oba; /* Band being optimised */ + int i, m, k; + + for (k = 0; k < (p->n * p->cord); k++) + dv[k] = 0.0; + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + int mo; + mppcol *c = &p->cols[i]; + double ov; + + /* Compute the tranfer corrected device values */ + /* and del in these values due to del in input parameters */ + for (m = 0; m < p->n; m++) { + tcnv[m] = icxdpTransFunc(&pv[m * p->cord], dtcnv_dpv[m], p->cord, c->nv[m]); + tcnv1[m] = 1.0 - tcnv[m]; + } + +#ifdef NEVER +{ +for (m = 0; m < p->n; m++) { + int ii; + double ttt; + + for (ii = 0; ii < p->cord; ii++) { + pv[m * p->cord + ii] += 1e-6; + ttt = (icxTransFunc(&pv[m * p->cord], p->cord, c->nv[m]) - tcnv[m])/1e-6; + pv[m * p->cord + ii] -= 1e-6; + printf("~1 chan %d cord %d is %f ref %f\n",m,ii,dtcnv_dpv[m][ii],ttt); + } +} +} +#endif + for (m = 0; m < p->n; m++) + ww[m] = 0.0; + + /* Lookup the shape values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + for (m = 0; m < p->n; m++) { + ww[m] += p->shape[m][k & ~(1<<m)][j] * vv; + /* Apply weighting to shape vertex value */ + } + } + + /* Compute del ww[m][m] for del tcnv[m] */ + for (m = 0; m < p->n; m++) { /* For each input channel were doing del of */ + for (mo = 0; mo < p->n; mo++) + dww_dtcnv[mo][m] = 0.0; + for (k = 0; k < p->nn; k++) { + int mm; + double vv = 1.0; + for (mm = 0; mm < p->n; mm++) { /* Compute weight for node k */ + if (m == mm) + continue; + if (k & (1 << mm)) + vv *= tcnv[mm]; + else + vv *= tcnv1[mm]; + } + for (mo = 0; mo < p->n; mo++) { /* For each output channel */ + double vvv = vv * p->shape[mo][k & ~(1<<mo)][j]; + if (k & (1 << m)) + dww_dtcnv[mo][m] += vvv; + else + dww_dtcnv[mo][m] -= vvv; + } + } + } + +#ifdef NEVER +{ + int mm; + double tww[MPP_MXINKS][MPP_MXINKS]; + + for (mm = 0; mm < p->n; mm++) { /* For del in each input value */ + tcnv[mm] += 1e-6; + + for (m = 0; m < p->n; m++) + tww[m][mm] = 0.0; + + /* Lookup the shape values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= (1.0 - tcnv[m]); + } + for (m = 0; m < p->n; m++) { + tww[m][mm] += p->shape[m][k & ~(1<<m)][j] * vv; + /* Apply weighting to shape vertex value */ + } + } + + tcnv[mm] -= 1e-6; + for (m = 0; m < p->n; m++) { + tww[m][mm] = (tww[m][mm] - ww[m])/1e-6; + } + } + + for (mm = 0; mm < p->n; mm++) { + for (m = 0; m < p->n; m++) { + printf("~1 ww[%d][%d] %f should be %f\n",mm, m, dww_dtcnv[mm][m],tww[mm][m]); + } + } +} +#endif /* NEVER */ + +#ifdef NEVER +{ + double itc[MPP_MXINKS]; + double ttc[MPP_MXINKS]; + double tww[MPP_MXINKS]; + + for (m = 0; m < p->n; m++) { + itc[m] = tcnv[m]; + } +#endif /* NEVER */ + + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double sv, dsv, vv = tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + tt = gg - gg * vv + 1.0; + sv = vv/tt; + dsv = (gg + 1.0)/(tt * tt); + } else { + tt = 1.0 - gg * vv; + sv = (vv - gg * vv)/tt; + dsv = (1.0 - gg)/(tt * tt); + } + tcnv[m] = sv; + tcnv1[m] = 1.0 - sv; + dtcnv_tc[m] = dsv; /* del in tcnv[m] due to del in tcnv[m] */ + dtcnv_ww[m] = (vv * vv - vv)/(tt * tt); /* del in tcnv[m] due to del in ww[m] */ + } + + +#ifdef NEVER + /* Check derivatives */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double vv = itc[m]; /* Input value to be tweaked */ + + vv += 1e-6; + if (gg >= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + ttc[m] = (vv - tcnv[m])/1e-6; + } + + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double vv = itc[m]; /* Input value to be tweaked */ + + gg += 1e-6; + if (gg >= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + tww[m] = (vv - tcnv[m])/1e-6; + } + + + for (m = 0; m < p->n; m++) { + printf("~1 dtcnv_tc[%d] is %f should be %f\n",m, dtcnv_tc[m], ttc[m]); + } + for (m = 0; m < p->n; m++) { + printf("~1 dtcnv_ww[%d] is %f should be %f\n",m, dtcnv_ww[m], tww[m]); + } +} +#endif /* NEVER */ + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double vv = p->pc[k][j]; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + /* Compute del ov for del tcnv[m] */ + for (m = 0; m < p->n; m++) { + for (dov[m] = 0.0, k = 0; k < p->nn; k++) { + int mm; + double vv = p->pc[k][j]; + for (mm = 0; mm < p->n; mm++) { + if (m == mm) + continue; + if (k & (1 << mm)) + vv *= tcnv[mm]; + else + vv *= tcnv1[mm]; + } + if (k & (1 << m)) + dov[m] += vv; + else + dov[m] -= vv; + } + } + +#ifdef NEVER +{ + int mm; + double tdov[MPP_MXINKS]; + + for (mm = 0; mm < p->n; mm++) { + tcnv[mm] += 1e-6; + + for (tdov[mm] = 0.0, k = 0; k < p->nn; k++) { + double vv = p->pc[k][j]; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= (1.0 - tcnv[m]); + } + tdov[mm] += vv; + } + tcnv[mm] -= 1e-6; + tdov[mm] = (tdov[mm] - ov)/1e-6; + } + + for (m = 0; m < p->n; m++) { + printf("~1 dov[%d] is %f should be %f\n",m, tdov[m], dov[m]); + } +} +#endif /* NEVER */ + + /* Compute the band value error */ + ddov = dDE(ov); + ov = lDE(ov) - c->lband[j]; + + ddov *= 2.0 * ov; /* del in final ov due to del in raw ov */ + rv += ov * ov; + + /* Accumulate delta from input params to output */ + for (m = 0; m < p->n; m++) { + for (k = 0; k < p->cord; k++) { + double ttt; + int mm; + + /* delta via dww */ + for (ttt = 0.0, mm = 0; mm < p->n; mm++) + ttt += dov[mm] * dtcnv_ww[mm] * dww_dtcnv[mm][m] * dtcnv_dpv[m][k]; + + /* delta direct */ + dv[m * p->cord + k] += ddov * (ttt + dov[m] * dtcnv_tc[m] * dtcnv_dpv[m][k]); + } + } + } + + rv /= (double)p->nodp; + for (k = 0; k < (p->n * p->cord); k++) + dv[k] /= (double)p->nodp; + + /* Compute weighted magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + tt = 2.0/(double)(p->n); /* Common factor */ + for (smv = 0.0, m = 0; m < p->n; m++) { + double w; + for (k = 0; k < p->cord; k++) { + i = m * p->cord + k; + w = (k < 2) ? TRANS_BASE : k * TRANS_HBASE; /* Increase weight with harmonics */ + dv[i] += w * tt * pv[i]; /* Del in rv due to del in pv */ + smv += w * pv[i] * pv[i]; + } + } + smv /= (double)(p->n); + rv += smv; + +#ifdef DEBUG + printf("dfunc2 itt %d/%d band %d returning %f\n",p->oit,p->ott,j,rv); +#endif + return rv; +} + +#ifdef TESTDFUNC +/* Check that dfunc2 returns the right values */ +static void test_dfunc2( +mpp *p, +int nparms, +double *pv +) { + int i, j, k; + double refov; + double refdv[MPP_MXINKS * MPP_MXTCORD]; + double ov; + double dv[MPP_MXINKS * MPP_MXTCORD]; + + /* Create reference dvs */ + refov = efunc2((void *)p, pv); + for (i = 0; i < nparms; i++) { + pv[i] += 1e-6; + refdv[i] = (efunc2((void *)p, pv) - refov)/1e-6; + pv[i] -= 1e-6; + } + + /* Check dfunc2 */ + ov = dfunc2((void *)p, dv, pv); + + printf("~#############################################\n"); + printf("~! Check dfunc2, ov %f, refov %f\n",ov, refov); + for (i = 0; i < nparms; i++) { + printf("~1 Parm %d val %f, ref %f\n",i,dv[i],refdv[i]); + } +} +#endif /* TESTDFUNC */ + + +/* Setup test point data ready for efunc3 on the given oba */ +static void sfunc3(mpp *p) { + double pcnv[MPP_MXCCOMB]; /* Interpolation combination values */ + int j = p->oba; /* Band being optimised */ + int i, k, m; + + /* Setup the correct per patch value fcnv values */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + + /* Compute the tranfer corrected device values */ + for (m = 0; m < p->n; m++) + c->tcnv[m] = icxTransFunc(p->tc[m][j], p->cord, c->nv[m]); + + /* Compute combination values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= c->tcnv[m]; + else + vv *= (1.0 - c->tcnv[m]); + } + pcnv[k] = vv; + } + + /* Compute shape weighting values */ + for (k = 0; k < p->nnn2; k++) { + int m = p->c2f[k].ink; + int n = p->c2f[k].comb; + /* Compress full interpolation one dimension to */ + /* exclude ink value being tweaked by shape */ + c->fcnv[k] = pcnv[n & ~(1<<m)] + pcnv[n | (1<<m)]; + } + } +} + + +/* Optimise all shape parameters simultaniously to minimise a particular bands error */ +/* Assume test point tcnv and pcnv are setup for pre-shape values */ +static double efunc3(void *adata, double pv[]) { + mpp *p = (mpp *)adata; + double smv, rv = 0.0; + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int j = p->oba; /* Band being optimised */ + int n1 = p->n - 1; + int i, m, k; + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double ov; + + for (m = 0; m < p->n; m++) + ww[m] = 0.0; + + /* Interpolate the per ink shape values for this set of input values */ + for (k = 0; k < p->nnn2; k++) { /* For each ink & interp vertex */ + m = k >> n1; /* Corresponding ink channel */ + ww[m] += pv[k] * c->fcnv[k]; /* Apply weighting to shape vertex value */ + } + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double vv = c->tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + tcnv[m] = vv; + tcnv1[m] = 1.0 - vv; + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double vv = p->pc[k][j]; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + /* Compute the band value error */ + ov = lDE(ov) - c->lband[j]; + rv += ov * ov; + } + + rv /= (double)p->nodp; + + /* Compute average magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + for (smv = 0.0, m = 0; m < p->nnn2; m++) { + smv += pv[m] * pv[m]; + } + smv /= (double)(p->nnn2); + rv += SHAPE_PMW * smv; + +#ifdef DEBUG + printf("efunc3 itt %d/%d band %d (smv %f) returning %f\n",p->oit,p->ott,j,smv,rv); +#endif + return rv; +} + +/* Return the gradient of the minimisation function at the given location, */ +/* as well as the function value at this location. */ +static double dfunc3(void *adata, double dv[], double pv[]) { + mpp *p = (mpp *)adata; + double smv, tt, rv = 0.0; + double dtcnv[MPP_MXINKS]; /* Derivative of transfer curve corrected device values */ + double dov[MPP_MXINKS]; /* Derivative of output interpolation device values */ + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int j = p->oba; /* Band being optimised */ + int n1 = p->n - 1; + int i, m, k; + + for (k = 0; k < p->nnn2; k++) + dv[k] = 0.0; /* Start with 0 delta */ + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double ov, ddov; + + for (m = 0; m < p->n; m++) + ww[m] = 0.0; + + /* Interpolate the per ink shape values for this set of input values */ + for (k = 0; k < p->nnn2; k++) { + m = k >> n1; /* Corresponding ink channel */ + ww[m] += pv[k] * c->fcnv[k]; /* Apply weighting to shape vertex value */ + } + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double sv, vv = c->tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + tt = gg - gg * vv + 1.0; + sv = vv/tt; + } else { + tt = 1.0 - gg * vv; + sv = (vv - gg * vv)/tt; + } + tcnv[m] = sv; + tcnv1[m] = 1.0 - sv; + dtcnv[m] = (vv * vv - vv)/(tt * tt); /* del in tcnv[m] due to del in ww[m] */ + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double vv = p->pc[k][j]; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + /* Compute del ov[m] for del tcnv[m] */ + for (m = 0; m < p->n; m++) { + for (dov[m] = 0.0, k = 0; k < p->nn; k++) { + int mm; + double vv = p->pc[k][j]; + for (mm = 0; mm < p->n; mm++) { + if (m == mm) + continue; + if (k & (1 << mm)) + vv *= tcnv[mm]; + else + vv *= tcnv1[mm]; + } + if (k & (1 << m)) + dov[m] += vv; + else + dov[m] -= vv; + } + dov[m] *= dtcnv[m]; /* del ov[m] due to del in ww[m] */ + } + + /* Compute the band value error */ + ddov = dDE(ov); + ov = lDE(ov) - c->lband[j]; + + ddov *= 2.0 * ov; /* del in final ov due to del in raw ov */ + rv += ov * ov; + + for (k = 0; k < p->nnn2; k++) { + m = k >> n1; /* Corresponding ink channel */ + dv[k] += ddov * dov[m] * c->fcnv[k]; + } + } + + rv /= ((double)p->nodp); + + for (k = 0; k < p->nnn2; k++) + dv[k] /= ((double)p->nodp); + + /* Compute average magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + tt = SHAPE_PMW * 2.0/(double)(p->nnn2); /* Common factor */ + for (smv = 0.0, k = 0; k < p->nnn2; k++) { + dv[k] += tt * pv[k]; /* del rv due to del pv[k] */ + smv += pv[k] * pv[k]; + } + smv /= (double)(p->nnn2); + rv += SHAPE_PMW * smv; + +#ifdef DEBUG + printf("dfunc3 itt %d/%d band %d (smv %f) returning %f\n",p->oit,p->ott,j,smv,rv); +#endif + return rv; +} + +#ifdef TESTDFUNC +/* Check that dfunc3 returns the right values */ +static void test_dfunc3( +mpp *p, +int nparms, +double *pv +) { + int i, j, k; + double refov; + double refdv[MPP_MXINKS * MPP_MXCCOMB/2]; + double ov; + double dv[MPP_MXINKS * MPP_MXCCOMB/2]; + + /* Create reference dvs */ + refov = efunc3((void *)p, pv); + for (k = 0; k < nparms; k++) { + pv[k] += 1e-6; + dv[k] = (efunc3((void *)p, pv) - refov)/1e-6; + pv[k] -= 1e-6; + } + + /* Check dfunc3 */ + ov = dfunc3((void *)p, dv, pv); + + printf("~#############################################\n"); + printf("~! Check dfunc3, ov %f, refov %f\n",ov, refov); + for (i = 0; i < nparms; i++) { + printf("~1 Parm %d val %f, ref %f\n",i,dv[i],refdv[i]); + } +} +#endif /* TESTDFUNC */ + + +/* Setup test point data ready for efunc4 on the given oba */ +static void sfunc4(mpp *p) { + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int j = p->oba; /* Band being optimised */ + int i, k, m; + + /* Setup the correct per patch value fcnv values */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + + /* Compute the tranfer corrected device values */ + for (m = 0; m < p->n; m++) { + tcnv[m] = icxTransFunc(p->tc[m][j], p->cord, c->nv[m]); + tcnv1[m] = 1.0 - tcnv[m]; + } + + for (m = 0; m < p->n; m++) + ww[m] = 0.0; + + /* Lookup the interpolated shape values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + for (m = 0; m < p->n; m++) { + ww[m] += p->shape[m][k & ~(1<<m)][j] * vv; + /* Apply weighting to shape vertex value */ + } + } + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double vv = tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + tcnv[m] = vv; + tcnv1[m] = 1.0 - vv; + } + + /* Compute the shape adjusted primary combination values */ + for (k = 0; k < p->nn; k++) { + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + c->pcnv[k] = vv; + } + } +} + +/* Optimise all vertex values simultaniously to minimise a particular bands error */ +/* Assume test point tcnv and pcnv are setup for post-shape values */ +static double efunc4(void *adata, double pv[]) { + mpp *p = (mpp *)adata; + double smv = 0.0, rv = 0.0; + int j = p->oba; /* Band being optimised */ + int i, k; + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double ov = 0.0; + + for (k = 0; k < p->nn; k++) { + double pp = pv[k]; + + if (pp < 0.0) /* Stop non-real values */ + rv += -5000.0 * pp; + ov += c->pcnv[k] * pp; + } + + /* Compute the band value error */ + ov = lDE(ov) - c->lband[j]; + rv += ov * ov; + } + + rv /= p->nodp; + + /* Compute anchor point error */ + for (smv = 0.0, i = 0; i < p->nn; i++) { + double tt = lDE(pv[i]) - p->lpca[i][j]; + smv += tt * tt; + } + smv /= (double)p->nn; + +//printf("~1 anchor error = %f\n",smv); + rv += COMB_PMW * smv; + +#ifdef DEBUG + printf("efunc4 itt %d/%d band %d returning %f\n",p->oit,p->ott,j,rv); +#endif + return rv; +} + +/* Return the gradient of the minimisation function at the given location, */ +/* as well as the function value at this location. */ +static double dfunc4(void *adata, double dv[], double pv[]) { + mpp *p = (mpp *)adata; + double smv = 0.0, rv = 0.0; + double drv[MPP_MXCCOMB]; /* Delta in rv */ + int j = p->oba; /* Band being optimised */ + int i, k; + + for (k = 0; k < p->nn; k++) { + dv[k] = 0.0; + drv[k] = 0.0; + } + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double ov = 0.0, ddov; + + for (k = 0; k < p->nn; k++) { + double pp = pv[k]; + + if (pp < 0.0) { /* Stop non-real values */ + rv += -5000.0 * pp; + drv[k] += -5000.0; + } + ov += c->pcnv[k] * pp; + } + + /* Compute the band value error */ + ddov = dDE(ov); + ov = lDE(ov) - c->lband[j]; + + ddov *= 2.0 * ov; /* del in final ov due to del in raw ov */ + rv += ov * ov; + + for (k = 0; k < p->nn; k++) { + dv[k] += ddov * c->pcnv[k]; + } + } + + rv /= p->nodp; + for (k = 0; k < p->nn; k++) + dv[k] /= ((double)p->nodp); + + /* Compute anchor point error */ + for (smv = 0.0, k = 0; k < p->nn; k++) { + double tt, ddov; + ddov = dDE(pv[k]); + tt = lDE(pv[k]) - p->lpca[k][j]; + ddov *= COMB_PMW * 2.0 * tt/(double)p->nn; + dv[k] += ddov; + smv += tt * tt; + } + smv /= (double)p->nn; + +//printf("~1 anchor error = %f\n",smv); + rv += COMB_PMW * smv; + + for (k = 0; k < p->nn; k++) /* Add any < 0 contribution */ + dv[k] += drv[k]; + +#ifdef DEBUG + printf("dfunc4 itt %d/%d band %d returning %f\n",p->oit,p->ott,j,rv); +#endif + return rv; +} + +#ifdef TESTDFUNC +/* Check that dfunc4 returns the right values */ +static void test_dfunc4( +mpp *p, +int nparms, +double *pv +) { + int i, j, k; + double refov; + double refdv[MPP_MXCCOMB]; + double ov; + double dv[MPP_MXCCOMB]; + + /* Create reference dvs */ + refov = efunc4((void *)p, pv); + for (i = 0; i < nparms; i++) { + pv[i] += 1e-9; + refdv[i] = (efunc4((void *)p, pv) - refov)/1e-9; + pv[i] -= 1e-9; + } + + /* Check dfunc4 */ + ov = dfunc4((void *)p, dv, pv); + + printf("~#############################################\n"); + printf("~! Check dfunc4, ov %f, refov %f\n",ov, refov); + for (i = 0; i < nparms; i++) { + printf("~1 Parm %d val %f, ref %f\n",i,dv[i],refdv[i]); + } +} +#endif /* TESTDFUNC */ + +/* ---------------------------------------------------- */ +/* Optimise whole model in one go */ + +#ifdef BIGBANG + +/* Optimise all parameters simultaniously to minimise a particular bands error */ +static double efunc0(void *adata, double pv[]) { + mpp *p = (mpp *)adata; + double *pv2, *pv3, *pv4; /* Pointers to each group of parameters */ + double smv, rv = 0.0; + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int j = p->oba; /* Band being optimised */ + int i, m, k; + + pv2 = pv; + pv3 = pv + (p->n * p->cord); + pv4 = pv + (p->n * p->cord) + p->nnn2; + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double ov; + + /* Compute the tranfer corrected device values */ + for (m = 0; m < p->n; m++) { + tcnv[m] = icxTransFunc(&pv2[m * p->cord], p->cord, c->nv[m]); + tcnv1[m] = 1.0 - tcnv[m]; + } + + for (m = 0; m < p->n; m++) + ww[m] = 0.0; + + /* Lookup the shape values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + for (m = 0; m < p->n; m++) { + ww[m] += pv3[p->f2c[m][k & ~(1<<m)]] * vv; + /* Apply weighting to shape vertex value */ + } + } + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double vv = tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + tcnv[m] = vv; + tcnv1[m] = 1.0 - vv; + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double vv = pv4[k]; + if (vv < 0.0) /* Stop non-real values */ + rv += -5000.0 * vv; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + /* Compute the band value error */ + ov = lDE(ov) - c->lband[j]; + rv += ov * ov; + } + + rv /= (double)p->nodp; + + /* Compute weighted magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + for (smv = 0.0, m = 0; m < p->n; m++) { + double w; + for (k = 0; k < p->cord; k++) { + i = m * p->cord + k; + w = (k < 2) ? TRANS_BASE : k * TRANS_HBASE; /* Increase weight with harmonics */ + smv += w * pv2[i] * pv2[i]; + } + } + smv /= (double)(p->n); + rv += smv; + + /* Compute average magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + for (smv = 0.0, m = 0; m < p->nnn2; m++) { + smv += pv3[m] * pv3[m]; + } + smv /= (double)(p->nnn2); + rv += SHAPE_PMW * smv; + + /* Compute anchor point error */ + for (smv = 0.0, i = 0; i < p->nn; i++) { + double tt = lDE(pv4[i]) - p->lpca[i][j]; + smv += tt * tt; + } + smv /= (double)p->nn; + rv += COMB_PMW * smv; + +#ifdef DEBUG + printf("efunc0 itt %d/%d band %d returning %f\n",p->oit,p->ott,j,rv); +#endif + return rv; +} + +/* Return the gradient of the minimisation function at the given location, */ +/* as well as the function value at this location. */ +static double dfunc0(void *adata, double dv[], double pv[]) { + mpp *p = (mpp *)adata; + double *pv2, *pv3, *pv4; /* Pointers to each group of parameters */ + double *dv2, *dv3, *dv4; /* Pointers to each group of derivatives */ + double smv, tt, rv = 0.0; + + double dtcnv_dpv2[MPP_MXINKS][MPP_MXTCORD]; /* Del in tcnv[m] due to del in pv2 */ + double dww_dpv3[MPP_MXINKS][MPP_MXINKS * MPP_MXCCOMB/2]; /* Del in ww[m] due to del in pv3 */ + double dov_dpv4[MPP_MXCCOMB]; /* Delta in ov due to delta in pv4 */ + double drv4[MPP_MXCCOMB]; /* Delta in rv due to error in pv4 */ + + double dww_dtcnv[MPP_MXINKS][MPP_MXINKS]; /* Del in ww[m] due to del in tcnv[m] */ + double dtcnv_tc[MPP_MXINKS]; /* Del in tcnv'[m] due to del in tcnv[m] */ + double dtcnv_ww[MPP_MXINKS]; /* Del in tcnv'[m] due to del in ww[m] */ + double dov_dcnv[MPP_MXINKS]; /* Del of ov due to del in tcnv'[m] */ + double ddov; /* Del in final ov due to del in raw ov */ + + double tcnv[MPP_MXINKS]; /* Transfer curve corrected device values */ + double tcnv1[MPP_MXINKS]; /* 1.0 - Transfer curve corrected device values */ + double ww[MPP_MXINKS]; /* Interpolated tweak params for each channel */ + int j = p->oba; /* Band being optimised */ + int i, m, k; + + pv2 = pv; + pv3 = pv + (p->n * p->cord); + pv4 = pv + (p->n * p->cord) + p->nnn2; + dv2 = dv; + dv3 = dv + (p->n * p->cord); + dv4 = dv + (p->n * p->cord) + p->nnn2; + + for (k = 0; k < (p->n * p->cord); k++) + dv2[k] = 0.0; + for (k = 0; k < p->nnn2; k++) + dv3[k] = 0.0; + for (k = 0; k < p->nn; k++) + dv4[k] += 0.0; + + /* For each test point */ + for (i = 0; i < p->nodp; i++) { + int mo; + mppcol *c = &p->cols[i]; + double ov; + + + /* Compute the tranfer corrected device values */ + /* and del in these values due to del in input parameters */ + for (m = 0; m < p->n; m++) { + tcnv[m] = icxdpTransFunc(&pv2[m * p->cord], dtcnv_dpv2[m], p->cord, c->nv[m]); + tcnv1[m] = 1.0 - tcnv[m]; + } + + for (m = 0; m < p->n; m++) { + ww[m] = 0.0; + for (k = 0; k < p->nnn2; k++) + dww_dpv3[m][k] = 0.0; + } + + /* Lookup the shape values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + for (m = 0; m < p->n; m++) { + int pix = p->f2c[m][k & ~(1<<m)]; + ww[m] += pv3[pix] * vv; /* Apply weighting to shape vertex value */ + dww_dpv3[m][pix] += vv; /* Delta to ww[m] from del in pv3[pix] */ + } + } + + /* Compute del ww[m][m] for del tcnv[m] */ + for (m = 0; m < p->n; m++) { /* For each input channel were doing del of */ + for (mo = 0; mo < p->n; mo++) + dww_dtcnv[mo][m] = 0.0; + for (k = 0; k < p->nn; k++) { + int mm; + double vv = 1.0; + for (mm = 0; mm < p->n; mm++) { /* Compute weight for node k */ + if (m == mm) + continue; + if (k & (1 << mm)) + vv *= tcnv[mm]; + else + vv *= tcnv1[mm]; + } + for (mo = 0; mo < p->n; mo++) { /* For each output channel */ + double vvv = vv * pv3[p->f2c[mo][k & ~(1<<mo)]]; + if (k & (1 << m)) + dww_dtcnv[mo][m] += vvv; + else + dww_dtcnv[mo][m] -= vvv; + } + } + } + + /* Apply the shape values to adjust the primaries */ + for (m = 0; m < p->n; m++) { + double gg = ww[m]; /* Curve adjustment */ + double sv, dsv, vv = tcnv[m]; /* Input value to be tweaked */ + if (gg >= 0.0) { + tt = gg - gg * vv + 1.0; + sv = vv/tt; + dsv = (gg + 1.0)/(tt * tt); + } else { + tt = 1.0 - gg * vv; + sv = (vv - gg * vv)/tt; + dsv = (1.0 - gg)/(tt * tt); + } + tcnv[m] = sv; + tcnv1[m] = 1.0 - sv; + dtcnv_tc[m] = dsv; /* del in tcnv[m] due to del in tcnv[m] */ + dtcnv_ww[m] = (vv * vv - vv)/(tt * tt); /* del in tcnv[m] due to del in ww[m] */ + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double pp, vv; + for (vv = 1.0, m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + dov_dpv4[k] = vv; + + pp = pv4[k]; + if (pp < 0.0) { /* Stop non-real values */ + rv += -5000.0 * pp; + drv4[k] += -5000.0; + } + ov += vv * pp; + } + + /* Compute del ov for del tcnv[m] */ + for (m = 0; m < p->n; m++) { + for (dov_dcnv[m] = 0.0, k = 0; k < p->nn; k++) { + int mm; + double vv = p->pc[k][j]; + for (mm = 0; mm < p->n; mm++) { + if (m == mm) + continue; + if (k & (1 << mm)) + vv *= tcnv[mm]; + else + vv *= tcnv1[mm]; + } + if (k & (1 << m)) + dov_dcnv[m] += vv; + else + dov_dcnv[m] -= vv; + } + } + + /* Compute the band value error */ + ddov = dDE(ov); + ov = lDE(ov) - c->lband[j]; + + ddov *= 2.0 * ov; /* del in final ov due to del in raw ov */ + rv += ov * ov; + + /* Accumulate delta from input params to output */ + for (m = 0; m < p->n; m++) { + for (k = 0; k < p->cord; k++) { + double ttt; + int mm; + + /* delta via dww */ + for (ttt = 0.0, mm = 0; mm < p->n; mm++) + ttt += dov_dcnv[mm] * dtcnv_ww[mm] * dww_dtcnv[mm][m] * dtcnv_dpv2[m][k]; + + /* delta direct */ + dv2[m * p->cord + k] += ddov * (ttt + dov_dcnv[m] * dtcnv_tc[m] * dtcnv_dpv2[m][k]); + } + } + for (k = 0; k < p->nnn2; k++) { + double ttt; + + /* delta via dww */ + for (ttt = 0.0, m = 0; m < p->n; m++) + ttt += dov_dcnv[m] * dtcnv_ww[m] * dww_dpv3[m][k]; + dv3[k] += ddov * ttt; + } + for (k = 0; k < p->nn; k++) { + dv4[k] += ddov * dov_dpv4[k]; + } + } + + rv /= (double)p->nodp; + for (k = 0; k < (p->n * p->cord); k++) + dv2[k] /= (double)p->nodp; + for (k = 0; k < p->nnn2; k++) + dv3[k] /= ((double)p->nodp); + for (k = 0; k < p->nn; k++) + dv4[k] /= ((double)p->nodp); + + /* Compute weighted magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + tt = 2.0/(double)(p->n); /* Common factor */ + for (smv = 0.0, m = 0; m < p->n; m++) { + double w; + for (k = 0; k < p->cord; k++) { + i = m * p->cord + k; + w = (k < 2) ? TRANS_BASE : k * TRANS_HBASE; /* Increase weight with harmonics */ + dv2[i] += w * tt * pv2[i]; /* Del in rv due to del in pv */ + smv += w * pv2[i] * pv2[i]; + } + } + smv /= (double)(p->n); + rv += smv; + + /* Compute average magnitude of shaper parameters squared */ + /* to minimise unconstrained "wiggles" */ + tt = SHAPE_PMW * 2.0/(double)(p->nnn2); /* Common factor */ + for (smv = 0.0, k = 0; k < p->nnn2; k++) { + dv3[k] += tt * pv3[k]; /* del rv due to del pv[k] */ + smv += pv3[k] * pv3[k]; + } + smv /= (double)(p->nnn2); + rv += SHAPE_PMW * smv; + + /* Compute anchor point error */ + for (smv = 0.0, k = 0; k < p->nn; k++) { + double tt, ddov; + ddov = dDE(pv4[k]); + tt = lDE(pv4[k]) - p->lpca[k][j]; + ddov *= COMB_PMW * 2.0 * tt/(double)p->nn; + dv4[k] += ddov; + smv += tt * tt; + } + smv /= (double)p->nn; + rv += COMB_PMW * smv; + + for (k = 0; k < p->nn; k++) /* Add any < 0 contribution */ + dv4[k] += drv4[k]; + +#ifdef DEBUG + printf("dfunc0 itt %d/%d band %d returning %f\n",p->oit,p->ott,j,rv); +#endif + return rv; +} + +#ifdef TESTDFUNC +/* Check that dfunc0 returns the right values */ +static void test_dfunc0( +mpp *p, +int nparms, +double *pv +) { + int i, j, k; + double refov; + double refdv[MPP_MXPARMS]; + double ov; + double dv[MPP_MXPARMS]; + + /* Create reference dvs */ + refov = efunc0((void *)p, pv); + for (i = 0; i < nparms; i++) { + pv[i] += 1e-9; + refdv[i] = (efunc0((void *)p, pv) - refov)/1e-9; + pv[i] -= 1e-9; + } + + /* Check dfunc0 */ + ov = dfunc0((void *)p, dv, pv); + + printf("~#############################################\n"); + printf("~! Check dfunc0, ov %f, refov %f\n",ov, refov); + for (i = 0; i < nparms; i++) { + printf("~1 Parm %d val %f, ref %f\n",i,dv[i],refdv[i]); + } +} +#endif /* TESTDFUNC */ +#endif /* BIGBANG */ + +#ifdef ISHAPE +/* ----------------------------------------------------- */ +/* Create an initial solution for shape parameters */ +/* by solving least squares using SVD */ + +/* Optimise device values to match test points Lab */ +/* while minimising the difference between the device */ +/* values and the current model device values in tcnv */ +static double efuncS(void *adata, double pv[]) { + double tcnv[MPP_MXINKS]; /* Shape tweaked input values */ + double tcnv1[MPP_MXINKS]; + mpp *p = (mpp *)adata; + double smv, rv = 0.0; + mppcol *c = p->otp; + int j = p->oba; + int k, m; + + double ov; + + for (k = 0; k < p->n; k++) { + double gg = pv[k]; /* Curve adjustment */ + double vv = c->tcnv[k]; /* Input value to be tweaked */ + if (gg >= 0.0) { + vv = vv/(gg - gg * vv + 1.0); + } else { + vv = (vv - gg * vv)/(1.0 - gg * vv); + } + tcnv[k] = vv; + tcnv1[k] = 1.0 - vv; + } + + /* Compute the primary combination values */ + for (ov = 0.0, k = 0; k < p->nn; k++) { + double vv = p->pc[k][j]; + for (m = 0; m < p->n; m++) { + if (k & (1 << m)) + vv *= tcnv[m]; + else + vv *= tcnv1[m]; + } + ov += vv; + } + + /* Compute the band value error */ + ov = lDE(ov) - c->lband[j]; + rv += ov * ov; + + /* Compute magnitude of shape params */ + for (smv = 0.0, k = 0; k < p->n; k++) { + smv += pv[k] * pv[k]; + } + rv += SHAPE_PMW * 0.1 * smv/(double)p->n; /* Don't worry about this here - SVD will cope */ + +//printf("~1 efuncS %f %f %f %f returning %f\n",pv[0],pv[1],pv[2],pv[3],rv); + return rv; +} + + +static void ishape( +mpp *p, +int j /* Band being initialised */ +) { + int i, k, m; + double **a; /* A[0..m-1][0..n-1] input A[][], will return U[][] */ + double *b; /* B[0..m-1] Right hand side of equation, return solution */ + + p->oba = j; + + if (p->nodp < (p->nn/2)) /* Nicer to return error ?? */ + error ("Not enough data points to solve shaper parameters"); + +//printf("~1 ishape called for band %d\n",j); + + /* First step is to compute the ideal shape values */ + /* for each test point. */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double pv[MPP_MXINKS]; /* Parameter values */ + double sr[MPP_MXINKS]; /* search radius */ + + p->otp = c; + +//printf("~1 setting up point %d\n",i); + + /* Compute the tranfer corrected device values */ + for (m = 0; m < p->n; m++) + c->tcnv[m] = icxTransFunc(p->tc[m][j], p->cord, c->nv[m]); + + /* Compute pre-shape combination values */ + for (k = 0; k < p->nn; k++) { /* For each interp vertex */ + double vv; + for (vv = 1.0, m = 0; m < p->n; m++) { /* Compute weighting */ + if (k & (1 << m)) + vv *= c->tcnv[m]; + else + vv *= (1.0 - c->tcnv[m]); + } + c->pcnv[k] = vv; + } + + /* Do reverse lookup to find ideal post transfer device */ + /* shape values that are minimal, but will give the desired */ + /* band value */ + for (k = 0; k < p->n; k++) { + pv[k] = 0.0; + sr[k] = 0.5; + } + + if (powell(NULL, p->n, pv, sr, 0.001, 300, efuncS, (void *)p, NULL, NULL) != 0) { +//printf("~1 ishape Powell failed at point %d\n",i); + for (k = 0; k < p->n; k++) + pv[k] = 0.0; + } + +//printf("~1 got uncorrected device vals %f %f %f %f\n",c->nv[0],c->nv[1],c->nv[2],c->nv[3]); +//printf("~1 got shape target values %f %f %f %f\n",pv[0],pv[1],pv[2],pv[3]); + + /* put them leave them in scnv */ + for (k = 0; k < p->n; k++) { + c->scnv[k] = pv[k]; + } + } + + /* The second step is to use SVD to compute the per colorant */ + /* interpolation corner values that are the best least squares */ + /* solution to providing those shape parameters at each point. */ + + /* Allocate our SVD matricies */ + a = dmatrix(0, p->nodp-1 + p->nn/2, 0, p->nn/2-1); + b = dvector(0, p->nodp-1 + p->nn/2); + + for (m = 0; m < p->n; m++) { /* For each channel */ + int kk; + + /* Setup our matricies for the test point */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double ss; /* Weighting factor */ + + /* The theoretically correct weighting to get */ + /* better alighment with DEsq() measure is */ + /* ss = dDE(c->band[j]);, but this doesn't seem */ + /* to work well in practice. */ + + /* I don't quite understand why this weighting */ + /* seems to give a good result. */ + ss = c->tcnv[m] * (1.0 - c->tcnv[m]); /* Don't understand why this is so good */ + + for (kk = k = 0; k < p->nn; k++) { + if (k & (1<<m)) + continue; /* Invalid comb for this channel */ + + a[i][kk++] = ss * (c->pcnv[k] + c->pcnv[k | (1<<m)]); + } + b[i] = ss * c->scnv[m]; + } + + /* Add extra fake "data points" to weight the minimisation of */ + /* the shape parameter values */ + { + int ii; + double ss; + + /* We want the relative weight to be similar to that */ + /* aimed for in efunc3 */ + + // ~~~ 0.1 seems about right for 1150 ??? + // 1.0 for roland8 ?? +// ss = (double)p->nodp/((double)(p->nn/2)) * 1.0 * SHAPE_PMW; + ss = (double)p->nodp * 0.05 * SHAPE_PMW; + + for (ii = 0; ii < p->nn/2; ii++) { + for (kk = k = 0; k < p->nn; k++) { + if (k & (1<<m)) + continue; /* Invalid comb for this channel */ + + if (kk == ii) /* Weights only on diagonals */ + a[i+ii][kk] = ss; + else + a[i+ii][kk] = 0.0; + kk++; + } + b[i+ii] = 0.0; + } + } + + /* Solve the equation A.x = b using SVD */ + if (svdsolve(a, b, p->nodp + p->nn/2, p->nn/2) != 0) + error("ishape: SVD failed"); + + /* Copy the solution back */ + for (kk = k = 0; k < p->nn; k++) { + if (k & (1<<m)) + continue; /* Invalid comb for this channel */ + +//printf("~1 Calculated shape[%d][%d] = %f\n",m,k,b[kk]); + p->shape[m][k][j] = b[kk++]; + } + } + + free_dvector(b, 0, p->nodp-1 + p->nn/2); + free_dmatrix(a, 0, p->nodp-1 + p->nn/2, 0, p->nn/2-1); +//printf("~1 ishape is done for band %d\n",j); + +} +#endif /* ISHAPE */ + + +/* ----------------------------------------------------- */ +/* Routine to figure out a suitable black point */ + +/* Structure to hold optimisation information */ +typedef struct { + mpp *p; /* Lookup object */ + int n; /* Number of device channels */ + double ilimit; /* Ink limit */ + double p1[3]; /* white pivot point */ + double p2[3]; /* Point on vector towards black */ +} bfinds; + +/* Optimise device values to minimise L, while remaining */ +/* within the ink limit, and staying in line between p1 and p2 */ +static double efunc6(void *adata, double pv[]) { + bfinds *b = (bfinds *)adata; + double rv = 0.0; + double dv[MPP_MXINKS]; + double Lab[3]; + double lr, ta, tb, terr; /* L ratio, target a, target b, target error */ + double sum, ovr; + int j; + + /* See if over ink limit or outside device range */ + for (ovr = sum = 0.0, j = 0; j < b->n; j++) { + dv[j] = pv[j]; + if (dv[j] < 0.0) { + if (-dv[j] > ovr) + ovr = -dv[j]; + dv[j] = 0.0; + } else if (dv[j] > 1.0) { + if ((dv[j] - 1.0) > ovr) + ovr = dv[j] - 1.0; + dv[j] = 1.0; + } + sum += dv[j]; + } + + if (b->ilimit > 1e-4) { + sum -= b->ilimit; + if (sum < 0.0) + sum = 0.0; + } else { + sum = 0.0; + } + + /* Compute Lab value */ + forward(b->p, NULL, Lab, NULL, dv); + +#ifdef DEBUG + printf("p1 = %f %f %f, p2 = %f %f %f\n",b->p1[0],b->p1[1],b->p1[2],b->p2[0],b->p2[1],b->p2[2]); + printf("device value %f %f %f %f, Lab = %f %f %f\n",dv[0],dv[1],dv[2],dv[3],Lab[0],Lab[1],Lab[2]); +#endif + + /* Primary is to minimise L value */ + rv = Lab[0]; + + /* See how out of line from p1 to p2 we are */ + lr = (Lab[0] - b->p1[0])/(b->p2[0] - b->p1[0]); /* Distance towards p2 from p1 */ + ta = lr * (b->p2[1] - b->p1[1]) + b->p1[1]; /* Target a value */ + tb = lr * (b->p2[2] - b->p1[2]) + b->p1[2]; /* Target b value */ + + terr = (ta - Lab[1]) * (ta - Lab[1]) + + (tb - Lab[2]) * (tb - Lab[2]); + +#ifdef DEBUG + printf("target error %f\n",terr); +#endif + rv += 100.0 * terr; + +#ifdef DEBUG + printf("out of range error %f\n",ovr); + printf("over limit error %f\n",sum); +#endif + rv += 200 * (ovr + sum); + +#ifdef DEBUG + printf("black find tc ret %f\n",rv); +#endif + return rv; +} + +static void compute_wb( +mpp *p +) { + int j; + + /* Figure out the white and black points */ + if (p->imask & ICX_ADDITIVE) { /* Assume additive doesn't have an ink limit */ + + /* Simply lookup values from min and max device. */ + /* If we have ICX_WHITE, should this be treated the same as */ + /* subtractive black ???? */ + for (j = 0; j < p->n; j++) + p->white.nv[j] = 1.0; + forward(p, &p->white.band[3], NULL, p->white.band, p->white.nv); + + for (j = 0; j < p->n; j++) + p->black.nv[j] = 0.0; + forward(p, &p->black.band[3], NULL, p->black.band, p->black.nv); + + for (j = 0; j < p->n; j++) + p->kblack.nv[j] = 0.0; + forward(p, &p->kblack.band[3], NULL, p->kblack.band, p->kblack.nv); + + } else { /* Subtractive */ + bfinds bfs; + double sr[MPP_MXINKS]; /* search radius */ + double tt[MPP_MXINKS]; /* temporary */ + int trial; + double rv, brv; + int kbset = 0; + + /* Lookup white directly */ + for (j = 0; j < p->n; j++) + p->white.nv[j] = 0.0; + forward(p, &p->white.band[3], bfs.p1, p->white.band, p->white.nv); + + /* Choose a black direction */ + if (p->imask & ICX_BLACK) { + int bix; + bix = icx_ink2index(p->imask, ICX_BLACK); + + for (j = 0; j < p->n; j++) { + if (j == bix) + p->kblack.nv[j] = 1.0; + else + p->kblack.nv[j] = 0.0; + } + forward(p, &p->kblack.band[3], bfs.p2, p->kblack.band, p->kblack.nv); + kbset = 1; + + } else if ((p->imask & (ICX_CYAN | ICX_MAGENTA | ICX_YELLOW)) + == (ICX_CYAN | ICX_MAGENTA | ICX_YELLOW)) { + int cix, mix, yix; + cix = icx_ink2index(p->imask, ICX_CYAN); + mix = icx_ink2index(p->imask, ICX_MAGENTA); + yix = icx_ink2index(p->imask, ICX_YELLOW); + + for (j = 0; j < p->n; j++) { + if (j == cix) + p->kblack.nv[j] = 1.0; + else if (j == mix) + p->kblack.nv[j] = 1.0; + else if (j == yix) + p->kblack.nv[j] = 1.0; + else + p->kblack.nv[j] = 0.0; + } + forward(p, NULL, bfs.p2, NULL, p->kblack.nv); + + } else { /* Make direction parallel to L axis */ + bfs.p2[0] = 0.0; + bfs.p2[1] = bfs.p1[1]; + bfs.p2[2] = bfs.p1[2]; + } + bfs.p = p; + bfs.n = p->n; + bfs.ilimit = p->limit; + + /* Find the black point */ + /* Do several trials to avoid local minima */ + for (j = 0; j < p->n; j++) { + tt[j] = p->black.nv[j] = 0.5; /* Starting point */ + sr[j] = 0.1; + } + brv = 1e38; + for (trial = 0; trial < 20; trial++) { + + if (powell(&rv, p->n, tt, sr, 0.00001, 500, efunc6, (void *)&bfs, NULL, NULL) == 0) { + if (rv < brv) { + brv = rv; + for (j = 0; j < p->n; j++) + p->black.nv[j] = tt[j]; + } + } + for (j = 0; j < p->n; j++) { + tt[j] = p->black.nv[j] + d_rand(-0.3, 0.3); + if (tt[j] < 0.0) + tt[j] = 0.0; + else if (tt[j] > 1.0) + tt[j] = 1.0; + } + } + if (brv > 1000.0) + error ("mpp: Black point powell failed"); + + for (j = 0; j < p->n; j++) { /* Make sure device values are in range */ + if (p->black.nv[j] < 0.0) + p->black.nv[j] = 0.0; + else if (p->black.nv[j] > 1.0) + p->black.nv[j] = 1.0; + } + + /* Set black value */ + forward(p, &p->black.band[3], NULL, p->black.band, p->black.nv); + + /* Set K only if device doesn't have a K channel */ + if (kbset == 0) { + for (j = 0; j < p->n; j++) + p->kblack.nv[j] = p->black.nv[j]; + forward(p, &p->kblack.band[3], NULL, p->kblack.band, p->kblack.nv); + } + + if (p->verb) { + double Lab[3]; + + icmXYZ2Lab(&icmD50, Lab, p->white.band); + printf("White point %f %f %f [Lab %f %f %f]\n", + p->white.band[0],p->white.band[1],p->white.band[2], + Lab[0], Lab[1], Lab[2]); + + icmXYZ2Lab(&icmD50, Lab, p->black.band); + printf("Black point %f %f %f [Lab %f %f %f]\n", + p->black.band[0],p->black.band[1],p->black.band[2], + Lab[0], Lab[1], Lab[2]); + + icmXYZ2Lab(&icmD50, Lab, p->kblack.band); + printf("K only Black point %f %f %f [Lab %f %f %f]\n", + p->kblack.band[0],p->kblack.band[1],p->kblack.band[2], + Lab[0], Lab[1], Lab[2]); + } + } +} + +/* ===================================== */ + +/* Create the mpp from scattered data points */ +/* Returns nz on error */ +static int create( + mpp *p, /* This */ + int verb, /* Vebosity level, 0 = none */ + int quality, /* Profile quality, 0..3 */ + int display, /* non-zero if display device */ + double limit, /* Total ink limit (if not display) */ + inkmask devmask, /* Inkmask describing device colorspace */ + 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) */ + double norm, /* Normalising scale value for spectral values */ + instType itype, /* Spectral instrument type (if not display) */ + int nodp, /* Number of points */ + mppcol *points /* Array of input points */ +) { + int it, i, j, k; + double de, mxde; /* Average Delta E and maximum Delta E */ + double sde, mxsde; /* Average Spectral error and maximum spectral error */ + int mxtcord; /* maximum transfer curve order */ + int maxit; /* Maximum number of tuning itterations */ + double thr; /* Powell threshold multiplier at each tuning pass */ + int useshape; /* Make use of shaping parameters */ + int mode; /* Band scanning mode */ + + /* Convert quality into operation counts */ + switch (quality) { + case 0: /* Low */ + useshape = 0; + mxtcord = 3; + maxit = 2; + break; + case 1: + default: /* Medium */ + useshape = 1; + mxtcord = 4; + maxit = 2; + break; + case 2: /* High */ + useshape = 1; + mxtcord = 5; + maxit = 3; + break; + case 3: /* Ultra high */ + useshape = 1; + mxtcord = 10; /* Is more actually better ? */ + maxit = 8; + break; + case 99: /* Special, simple model */ + useshape = 0; + mxtcord = 1; + maxit = 4; + break; + } + + /* Setup the basic mpp information */ + p->verb = verb; + p->imask = devmask; + p->n = icx_noofinks(devmask); + p->nn = 1 << p->n; + p->nnn2 = p->n * p->nn/2; + + if (display) { + p->display = 1; + p->limit = p->n; + p->itype = instUnknown; + } else { + p->display = 0; + p->limit = limit; /* Record it here */ + p->itype = itype; + } + p->spec_n = spec_n; + p->spec_wl_short = spec_wl_short; + p->spec_wl_long = spec_wl_long; + p->nodp = nodp; + + /* MPP limit is less than XICC */ + if (p->n > MPP_MXINKS) { + p->errc = 1; + sprintf(p->err,"MPP Can't handle %d colorants",p->n); + return 1; + } + + /* MPP limit is less than XSPECT */ + if (spec_n > MPP_MXBANDS) { + p->errc = 1; + sprintf(p->err,"MPP Can't handle %d spectral bands",spec_n); + return 1; + } + + /* Take a copy of the data points */ + if ((p->cols = new_mppcols(p->nodp, p->n, p->spec_n)) == NULL) { + error("Malloc failed!"); + } + if ((new_mppcol(&p->white, p->n, p->spec_n)) != 0) { + error("Malloc failed!"); + } + if ((new_mppcol(&p->black, p->n, p->spec_n)) != 0) { + error("Malloc failed!"); + } + if ((new_mppcol(&p->kblack, p->n, p->spec_n)) != 0) { + error("Malloc failed!"); + } + + p->spmax = -1e6; + for (i = 0; i < p->nodp; i++) { + copy_mppcol(&p->cols[i], &points[i], p->n, p->spec_n); /* Copy structure */ + + /* Create Lab version */ + icmXYZ2Lab(&icmD50, p->cols[i].Lab, p->cols[i].band); + +#ifdef SHARPEN + XYZ2sharp(&p->cols[i].band[0], &p->cols[i].band[1], &p->cols[i].band[2]); +#endif + /* Normalise spectral values */ + for (j = 0; j < p->spec_n; j++) { + p->cols[i].band[3+j] /= norm; /* Normalise spectral value to range 0..1 */ + + if (p->cols[i].band[3+j] > p->spmax) /* Track maximum value */ + p->spmax = p->cols[i].band[3+j]; + } + + /* Compute L* type band target values */ + for (j = 0; j < (3+p->spec_n); j++) { + p->cols[i].lband[j] = lDE(p->cols[i].band[j]); + } + } + p->norm = 1.0; /* Internal norm is 1.0 */ + + /* Compute L* type band target values */ + for (i = 0; i < p->nodp; i++) { + for (j = 0; j < (3+p->spec_n); j++) { + p->cols[i].lband[j] = lDE(p->cols[i].band[j]); + } + } + /* Init transfer curve parameters of model */ + for (k = 0; k < p->n; k++) { /* For each ink */ + for (j = 0; j < (p->spec_n+3); j++) { /* For each band */ + for (i = 0; i < mxtcord; i++) { /* For each curve order */ + if (i == 0) + p->tc[k][j][i] = -1.6; /* Typical starting value */ + else + p->tc[k][j][i] = 0.0; /* Straight transfer curve */ + } + } + } + + + /* Initialise the primary colorant values */ + for (i = -1; i < p->n; i++) { + int ii, bk = 0; + double bdif = 1e6; + + if (i < 0) + ii = 0; + else + ii = 1 << i; + + /* Search the patch list to find the one closest to this colorant combination */ + for (k = 0; k < p->nodp; k++) { + double dif = 0.0; + for (j = 0; j < p->n; j++) { + double tt; + if (i == j) + tt = 1.0 - p->cols[k].nv[j]; + else + tt = p->cols[k].nv[j]; + dif += tt * tt; + } + if (dif < bdif) { /* best so far */ + bdif = dif; + bk = k; + if (dif < 0.001) + break; /* Don't bother looking further */ + } + } + + /* Put that sample patch in place as initial value */ + for (j = 0; j < (3+p->spec_n); j++) + p->pc[ii][j] = p->cols[bk].band[j]; +#ifdef DEBUG + printf("comb 0x%x XYZ is %f %f %f\n", ii, p->pc[ii][0], p->pc[ii][1], p->pc[ii][2]); +#endif + } + + /* Estimate primary combination values from primary values */ + { + double sm[3+MPP_MXBANDS]; /* Smallest reflection sample in this band */ + double ink[3+MPP_MXBANDS]; + + /* Search the patch list to find the smallest values for each band */ + /* we'll use this as a guide to the freznell reflection from the surface */ + for (j = 0; j < (3+p->spec_n); j++) { + + sm[j] = 1e6; + for (k = 0; k < p->nodp; k++) { + if (sm[j] > p->cols[k].band[j]) { + int m; + sm[j] = p->cols[k].band[j]; + ink[j] = 0.0; + for (m = 0; m < p->n; m++) + ink[j] += p->cols[k].nv[m]; + } + } + /* Adjust for error in freznell was estimated from */ + /* a low ink coverage */ + if (ink[j] < 4.0) { + sm[j] *= pow(0.775, 4.0 - ink[j]); + } +#ifdef DEBUG + printf("smallest value in band %d = %f from total ink %f\n",j,sm[j],ink[j]); +#endif + } + + for (i = 3; i < p->nn; i++) { + + if ((i & (i-1)) == 0) + continue; /* Skip primaries */ + for (j = 0; j < (3+p->spec_n); j++) { + int k; + double wh = p->pc[0][j]; /* white */ + double tr = 1.0; /* Trapping coefficient */ + if (wh < 0.01) + wh = 0.01; /* Guard against silliness */ + p->pc[i][j] = wh; /* Start with white */ + for (k = 0; k < p->n; k++) { + if (i & (1<<k)) { + double co = p->pc[1 << k][j]; /* colorant + paper reflectance */ + co = (co - sm[j])/wh; /* Colorant reflectance */ + co /= tr; /* Trapping reduction */ + if (co > 1.0) + co = 1.0; + else if (co < 0.0) + co = 0.0; + p->pc[i][j] *= co; /* Estimated combined reflectivity */ + + tr *= 0.90; /* Next inks effectiveness due to trapping */ + } + } + p->pc[i][j] = p->pc[i][j]; + p->pc[i][j] += sm[j]; /* Add freznell back in */ + } +#ifdef DEBUG + printf("comb 0x%x estimated XYZ = %f %f %f\n",i, p->pc[i][0], p->pc[i][1], p->pc[i][2]); +#endif + } + } + + /* Override the estimated primary combination values, if actual readings are available */ + for (i = 3; i < p->nn; i++) { + int k, bk = 0; + double bdif = 1e6; + + if ((i & (i-1)) == 0) + continue; /* Skip primaries */ + + /* Search the patch list to find the one closest to this colorant combination */ + for (k = 0; k < p->nodp; k++) { + double dif = 0.0; + for (j = 0; j < p->n; j++) { + double tt; + if (i & (1<<j)) + tt = 1.0 - p->cols[k].nv[j]; + else + tt = p->cols[k].nv[j]; + dif += tt * tt; + } + if (dif < bdif) { /* best so far */ + bdif = dif; + bk = k; + if (dif < 0.001) + break; /* Don't bother looking further */ + } + } + + if (bdif < 0.02) { + + /* Override the estimated combination values with the real values */ + for (j = 0; j < (3+p->spec_n); j++) { +#ifdef DEBUG + printf("comb 0x%x band %d was %f ",i, j, p->pc[i][j]); +#endif + p->pc[i][j] = p->cols[bk].band[j]; +#ifdef DEBUG + printf("best %d now %f\n",bk, p->pc[i][j]); +#endif + } + } + } + + /* These initial primary combination values become the anchors during optimisation */ + for (i = 0; i < p->nn; i++) { + for (j = 0; j < (3+p->spec_n); j++) + p->lpca[i][j] = lDE(p->pc[i][j]); + } + + /* Allocate and init shape related parameter space */ + init_shape(p); + + p->cord = 1; /* Start with only 1 order */ + +#ifndef DEBUG + if (p->verb) +#endif /* !DEBUG */ + { + deltae(p, &de, &mxde, &sde, &mxsde); + printf("Before optimising model have average dE of %f, max %f\n", de, mxde); + if (p->spec_n > 0) printf("and average spectral E of %f, max %f\n",sde, mxsde); + } + + /* - - - - - - - - - - - - - - - - - */ + +#ifndef NOPROCESS +#ifdef NEVER // Skip efunc1 passes for now. + /* Do initial fast pass of optimisations */ + /* using only first transfer curve order */ + for (it = 0, p->ott = 3; it < p->ott; it++) { + double resid; + + p->oit = it+1; + + /* First optimise each input channels transfer curve */ + for (k = 0; k < p->n; k++) { /* For each input channel */ + double tw = 0.0; /* Total weight */ + + p->och = k; /* device channel being optimised */ + + /* Setup the appropriate weights */ + for (i = 0; i < p->nodp; i++) { + mppcol *c = &p->cols[i]; + double ww; + int kk; + + for (ww = 0.0, kk = 0; kk < p->n; kk++) { + if (kk == k) + continue; /* Channel of interest can have any value */ + ww += c->nv[kk] * c->nv[kk]; + } + c->w = 1.0 - ww; + if (c->w < 0.0) + c->w = 0.0; + tw += c->w; + +//printf("~1 chan %d point %d weight %f\n",k,i,c->w); + } + + if (tw < 3.0) + error ("MPP - not enough weighted point"); + + for (j = 0; j < (3+p->spec_n); j++) { /* For all bands */ + double pv[MPP_MXTCORD]; /* Parameter values */ + double sr[MPP_MXTCORD]; /* search radius */ + + p->oba = j; /* Band being optimised */ + + if (p->verb) { + printf("Optimising device transfer curves channel %d band %d:\n",k,j); + } + + /* Get the current values */ + for (i = 0; i < p->cord; i++) { + pv[i] = p->tc[k][j][i]; + sr[i] = 0.05; + } + + sfunc1(p); /* Setup test point values for this chan and band */ + + if (powell(&resid, p->cord, pv, sr, 0.001, 100, efunc1, (void *)p, mppprog, (void *)p) != 0) + error ("Powell failed"); + + /* Put results back into place */ + for (i = 0; i < p->cord; i++) { + p->tc[k][j][i] = pv[i]; + } + } +#ifndef DEBUG + if (p->verb) +#endif /* !DEBUG */ + { + deltae(p, &de, &mxde, &sde, &mxsde); + printf("\nNow got avg dE of %f, max %f\n",de, mxde); + if (p->spec_n > 0) + printf("and avgerage spectral E of %f, max %f\n",sde, mxsde); + } + } + } +#endif /* NEVER */ + + p->cord = mxtcord; + + /* - - - - - - - - - - - - - - - - - */ + /* Fine tune all parameters in the model */ + for (mode = 0;;) { + double pv[MPP_MXPARMS]; /* Parameter values */ + double sr[MPP_MXPARMS]; /* search radius */ + int lj, yj = 0; /* Last band, peak Y band */ + + /* Decide which band to do next */ + if (mode == 0) { /* Start at the beginning */ + lj = -1; + if (p->spec_n == 0) { + mode = 3; /* Switch to doing Y values */ + j = 1; + } else { + /* Start at the peak Y bandwidth */ + j = (int)(p->spec_n * (555.0 - p->spec_wl_short) + /(p->spec_wl_long - p->spec_wl_short) + 0.5); + if (j < 0) + j = 0; + else if (j >= p->spec_n) + j = p->spec_n-1; + j += 3; + yj = j; + mode = 1; /* Switch to incrementing j */ + } + } else if (mode == 1) { /* Increment j */ + lj = j; + j++; + if (j >= (3+p->spec_n)) { /* We're finished moving up */ + lj = yj; + j = yj-1; + mode = 2; + if (j < 3) { /* Just in case */ + lj = yj; + j = 1; + mode = 3; + } + } + } else if (mode == 2) { /* Decrement j */ + lj = j; + j--; + if (j < 3) { + lj = yj; /* We know that spect is valid */ + j = 1; + mode = 3; /* Switch to Y mode */ + } + } else if (mode == 3) { /* Doing Y mode */ + if (spec_n > 0) { + /* Use spec at the peak X bandwidth */ + lj = (int)(p->spec_n * (600.0 - p->spec_wl_short) + /(p->spec_wl_long - p->spec_wl_short) + 0.5); + if (lj < 0) + lj = 0; + else if (lj >= p->spec_n) + lj = p->spec_n-1; + lj += 3; + } else { + lj = 1; + } + j = 0; /* Doing X mode */ + mode = 4; + } else if (mode == 4) { /* Doing X mode */ + if (spec_n > 0) { + /* Use spec at the peak Z bandwidth */ + lj = (int)(p->spec_n * (445.0 - p->spec_wl_short) + /(p->spec_wl_long - p->spec_wl_short) + 0.5); + if (lj < 0) + lj = 0; + else if (lj >= p->spec_n) + lj = p->spec_n-1; + lj += 3; + } else { + lj = 1; + } + j = 2; /* Doing Z mode */ + mode = 5; + } else { + break; /* we're now done */ + } + + p->oba = j; /* Band being optimised */ + + if (p->verb) + printf("Doing band %d, last band %d\n",j,lj); + + /* See if the last bands values are a good place to start */ + if (lj >= 0) { + double cval, pval, p0val; + + banderr(p, &cval, NULL, j); /* Current error */ + + /* Copy previous band transfer and shape into current band */ + for (k = 0; k < p->n; k++) + for (i = 0; i < p->cord; i++) { + pv[k * p->cord + i] = p->tc[k][j][i]; /* Save current for restore */ + p->tc[k][j][i] = p->tc[k][lj][i]; + } + for (i = 0; i < p->nnn2; i++) { + int m = p->c2f[i].ink; + int n = p->c2f[i].comb; + + sr[i] = p->shape[m][n][j]; + p->shape[m][n][j] = p->shape[m][n][lj]; + } + banderr(p, &pval, NULL, j); + + /* Try out the urrent order 0 transfer values with rest of transfer and shape */ + for (k = 0; k < p->n; k++) + p->tc[k][j][0] = pv[k * p->cord]; + banderr(p, &p0val, NULL, j); + + /* See which was best out of the three */ + if (pval >= cval && p0val >= pval) { /* Original was the best */ + for (k = 0; k < p->n; k++) + for (i = 0; i < p->cord; i++) + p->tc[k][j][i] = pv[k * p->cord + i]; /* Restore previous values */ + for (i = 0; i < p->nnn2; i++) { + int m = p->c2f[i].ink; + int n = p->c2f[i].comb; + p->shape[m][n][j] = sr[i]; /* Restore previous value */ + } +//printf("~1 Starting values were best (%f && %f > %f)\n",pval,p0val,cval); + } else if (p0val >= pval) { /* Copying all was best */ + + for (k = 0; k < p->n; k++) + p->tc[k][j][0] = p->tc[k][lj][0]; /* Back to order 0 values */ + +//printf("~1 copying all previous bands values was best (%f < %f && %f)\n",pval,cval,p0val); + } else { +//printf("~1 copying except order 0 values was best (%f < %f && %f)\n",p0val,cval,pval); + } + } + +#ifdef MULTIPASS /* Multipass in parts */ + for (it = 0, p->ott = maxit, thr = 1.0; it < maxit; it++, thr *= 0.2) { + double sde, mxsde; + double resid; + + p->oit = it+1; + + /* Optimise main transfer curve to minimise each bands error */ + /* Initially using only first transfer curve order */ + + if (p->verb) + printf("Fine tuning device transfer curves itteration %d\n",it); + + /* Get the current values */ + for (k = 0; k < p->n; k++) { + for (i = 0; i < p->cord; i++) { + pv[k * p->cord + i] = p->tc[k][j][i]; + sr[k * p->cord + i] = 0.05; + } + } +#ifdef TESTDFUNC + test_dfunc2(p, p->n * p->cord, pv); +#endif /* TESTDFUNC */ +#ifdef NODDV + if (powell(&resid, p->n * p->cord, pv, sr, thr * 0.01, 200, + efunc2, (void *)p, mppprog, (void *)p) != 0) + error ("Powell failed"); +#else /* !NODDV */ + if (conjgrad(&resid, p->n * p->cord, pv, sr, thr * 0.01, 200, + efunc2, dfunc2, (void *)p, mppprog, (void *)p)!= 0) + error ("ConjGrad failed"); +#endif /* !NODDV */ + + /* Put results back into place */ + for (k = 0; k < p->n; k++) { + for (i = 0; i < p->cord; i++) + p->tc[k][j][i] = pv[k * p->cord + i]; + } + +#ifndef DEBUG + if (p->verb) +#endif /* !DEBUG */ + { + banderr(p, &sde, &mxsde, j); /* Current error */ + printf("\nNow got avg E of %f, max %f for this band\n",sde, mxsde); + } + + p->cord = mxtcord; /* maximum transfer curve order after very first run */ + + /* Tune the shaping parameters */ + if (useshape) { + + if (p->verb) + printf("Tuning detailed shaping parameters itteration %d\n",it); + + sfunc3(p); /* Setup test point values for this band */ + + /* Get the current values */ + for (i = 0; i < p->nnn2; i++) { + int m = p->c2f[i].ink; + int n = p->c2f[i].comb; + + pv[i] = p->shape[m][n][j]; + sr[i] = 0.01; + } + +#ifdef TESTDFUNC + test_dfunc3(p, p->nnn2, pv); +#endif /* TESTDFUNC */ +#ifdef NODDV + if (powell(&resid, p->nnn2, pv, sr, thr * 0.05, 2000, + efunc3, (void *)p, mppprog, (void *)p) != 0) + error ("Powell failed"); + +#else /* !NODDV */ + if (conjgrad(&resid, p->nnn2, pv, sr, thr * 0.05, 2000, + efunc3, dfunc3, (void *)p, mppprog, (void *)p) != 0.0) + error ("ConjGrad failed"); +#endif /* !NODDV */ + + /* Put results back into place */ + for (i = 0; i < p->nnn2; i++) { + int m = p->c2f[i].ink; + int n = p->c2f[i].comb; + + p->shape[m][n][j] = pv[i]; +//printf("~1 shape[%d][%d] = %f\n",m,n,pv[i]); + } + +#ifndef DEBUG + if (p->verb) +#endif /* !DEBUG */ + { + banderr(p, &sde, &mxsde, j); /* Current error */ + printf("\nNow got avg E of %f, max %f for this band\n",sde, mxsde); + } +#ifdef DEBUG + dump_shape(p, 0, "After efunc3:"); +#endif /* DEBUG */ + + p->useshape = 1; /* Would be nice to flag this on a per band basis */ + } + + /* Tune the vertex parameters */ + + if (p->verb) + printf("Optimising device combination values itteration %d\n",it); + + sfunc4(p); /* Setup test point values for this band */ + + /* Get the current values */ + for (k = 0; k < p->nn; k++) { + pv[k] = p->pc[k][j]; + sr[k] = 0.01; + } + +#ifdef TESTDFUNC + test_dfunc4(p, p->nn, pv); +#endif /* TESTDFUNC */ +#ifdef NODDV + if (powell(&resid, p->nn, pv, sr, thr * 0.01, 500, + efunc4, (void *)p, mppprog, (void *)p) != 0) + error ("Powell failed"); +#else /* !NODDV */ + if (conjgrad(&resid, p->nn, pv, sr, thr * 0.01, 500, + efunc4, dfunc4, (void *)p, mppprog, (void *)p) != 0) + error ("ConjGrad failed"); +#endif /* !NODDV */ + + /* Put results back into place */ + for (k = 0; k < p->nn; k++) { + double pp = pv[k]; + if (pp < 0.0) + pp = 0.0; + p->pc[k][j] = pp; + } + +#ifndef DEBUG + if (p->verb) +#endif /* !DEBUG */ + { + banderr(p, &sde, &mxsde, j); /* Current error */ + printf("\nNow got avg E of %f, max %f for this band\n",sde, mxsde); + } + } +#endif /* MULTIPASS */ + +#ifdef BIGBANG + /* Do optimisation with one big bang */ + { + double resid; + double *pv2, *pv3, *pv4; /* Pointers to each group of parameters */ + double *sr2, *sr3, *sr4; /* Pointers to each group of search radius */ + int tparms = p->n * p->cord + p->nnn2 + p->nn; + + pv2 = pv; + pv3 = pv + (p->n * p->cord); + pv4 = pv + (p->n * p->cord) + p->nnn2; + sr2 = sr; + sr3 = sr + (p->n * p->cord); + sr4 = sr + (p->n * p->cord) + p->nnn2; + + /* Get the current transfer values */ + for (k = 0; k < p->n; k++) { + for (i = 0; i < p->cord; i++) { + pv2[k * p->cord + i] = p->tc[k][j][i]; + sr2[k * p->cord + i] = 0.005; + } + } + /* Get the current shaper values */ + for (i = 0; i < p->nnn2; i++) { + int m = p->c2f[i].ink; + int n = p->c2f[i].comb; + + pv3[i] = p->shape[m][n][j]; + sr3[i] = 0.005; + } + /* Get the current device combination values */ + for (k = 0; k < p->nn; k++) { + pv4[k] = p->pc[k][j]; + sr4[k] = 0.005; + } + +#ifdef TESTDFUNC + test_dfunc0(p, tparms, pv); +#endif /* TESTDFUNC */ + + if (conjgrad(&resid, tparms, pv, sr, 0.001, 4000, efunc0, dfunc0, (void *)p, + mppprog, (void *)p) != 0) + error ("ConjGrad failed"); + + /* Put results back into place */ + for (k = 0; k < p->n; k++) { + for (i = 0; i < p->cord; i++) + p->tc[k][j][i] = pv2[k * p->cord + i]; + } + for (i = 0; i < p->nnn2; i++) { + int m = p->c2f[i].ink; + int n = p->c2f[i].comb; + + p->shape[m][n][j] = pv3[i]; + } + for (k = 0; k < p->nn; k++) { + double pp = pv4[k]; + if (pp < 0.0) + pp = 0.0; + p->pc[k][j] = pp; + } + +#ifndef DEBUG + if (p->verb) +#endif /* !DEBUG */ + { + banderr(p, &sde, &mxsde, j); /* Current error */ + printf("\nNow got avg E of %f, max %f for this band\n",sde, mxsde); + } + } +#endif /* BIGBANG */ + } +#endif /* NOPROCESS */ + +#ifdef DEBUG + if (p->verb) +#endif /* !DEBUG */ + { + double de, mxde; + deltae(p, &de, &mxde, &sde, &mxsde); + printf("\nNow got avg dE of %f, max %f\n",de, mxde); + if (p->spec_n > 0) + printf("and avgerage spectral E of %f, max %f\n",sde, mxsde); + } + +#ifdef DOPLOT /* Plot the device curves */ + { +#define XRES 100 + double xx[XRES]; + double y1[XRES]; + + for (i = 0; i < (3+p->spec_n); i++) { + printf("Band %d:\n",i); + for (j = 0; j < p->n; j++) { + printf("Ink %d:\n",j); + + for (k = 0; k < XRES; k++) { + double x; + x = k/(double)(XRES-1); + xx[k] = x; + y1[k] = icxTransFunc(p->tc[j][i], p->cord, x); + } + do_plot(xx,y1,NULL,NULL,XRES); + } + } + } +#endif /* DOPLOT */ + + /* Figure out the white and black points */ + compute_wb(p); + + /* Done with our copy of the input points */ + free (p->cols); + p->nodp = 0; + p->cols = NULL; + + return 0; +} + + +#ifdef DEBUG + +/* Dump current shape params to a file */ +static void dump_shape(mpp *p, int first, char *title) { + int i,j,k; + FILE *df; + /* Some debug code */ + + if (first) { + if ((df = fopen("debug.txt","w")) == NULL) + error ("Failed to open debug.txt"); + } else { + if ((df = fopen("debug.txt","a")) == NULL) + error ("Failed to open debug.txt"); + } + + fprintf(df,"%s\n",title); + + for (i = 0; i < (3+p->spec_n); i++) { + fprintf(df,"Band %d:\n",i); + for (j = 0; j < p->n; j++) { + fprintf(df,"Ink %d:\n",j); + + for (k = 0; k < p->nn; k++) { + if ((k & (1<<j)) == 0) { + int m, n; + double val = 0.0; + fprintf(df,"Comb %d = %f\n", k, p->shape[j][k][i]); + } + } + fprintf(df,"\n"); + } + fprintf(df,"\n"); + } + fclose(df); +} + +/* Combine two first order shapers coefficients together */ +static double comb(double n1, double n2) { + double nn; + if (n1 > 0.0) + n1 = (n1+1.0); + else + n1 = 1.0/(1.0-n1); + + if (n2 > 0.0) + n2 = (n2+1.0); + else + n2 = 1.0/(1.0-n2); + + nn = n1 * n2; + + if (nn >= 1.0) + nn -= 1.0; + else + nn = 1.0-(1.0/nn); + + return nn; +} + +#endif /* DEBUG */ + + + + + + + + + + + + + + + + + + + + + + + + + |