summaryrefslogtreecommitdiff
path: root/profile/printcal.c
diff options
context:
space:
mode:
Diffstat (limited to 'profile/printcal.c')
-rw-r--r--profile/printcal.c2320
1 files changed, 2320 insertions, 0 deletions
diff --git a/profile/printcal.c b/profile/printcal.c
new file mode 100644
index 0000000..f9107eb
--- /dev/null
+++ b/profile/printcal.c
@@ -0,0 +1,2320 @@
+
+
+
+/*
+ * Argyll Color Correction System
+ * Print Device calibration curve generator.
+ *
+ * Author: Graeme W. Gill
+ * Date: 2008/3/3
+ *
+ * Copyright 1996-2008 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 program takes in the colorent wedge test chart
+ * points, and creates a set of per channel correction curves.
+ */
+
+/*
+ * TTBD:
+ * Allow auto max threshold to be scaled on command line ?
+ * ie. -m# set % to go below the default optimal maximum.
+ */
+
+
+/*
+ Additive spaces are handled by inverting the device values internally.
+ (Such a space should probably have ICX_INVERTED set as well, indicating
+ that the underlying device is actually subtractive.)
+ */
+
+#undef DEBUG
+
+#define verbo stdout
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <fcntl.h>
+#include <string.h>
+#include <time.h>
+#include "copyright.h"
+#include "aconfig.h"
+#include "cgats.h"
+#include "numlib.h"
+#include "sort.h"
+#include "rspl.h"
+#include "xicc.h"
+#include "plot.h"
+
+
+#define RSPLFLAGS (0 /* | RSPL_2PASSSMTH | RSPL_EXTRAFIT2 */)
+
+#define RSPLSMOOTH 2.0 /* RSPL Smoothness factor use on measured device points */
+
+#define TCURVESMOOTH 1.0 /* RSPL smoothness factor for target aim points */
+
+#define GRES 256 /* Rspl grid resolution */
+#define SLOPE_NORM 70.0 /* Normalized delta E for below thresholds */
+#define MIN_SLOPE_A 8.0 /* Criteria for Auto max, DE/dDev at max */
+#define MIN_SLOPE_O 3.0 /* Criteria for Auto max, min DE/dDev below max */
+
+#define CAL_RES 256 /* Resolution saved to .cal file */
+
+#define PRES 256 /* Plotting resolution */
+
+void usage(char *diag, ...) {
+ int i;
+ fprintf(stderr,"Create printer calibration, Version %s\n",ARGYLL_VERSION_STR);
+ fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
+ if (diag != NULL) {
+ va_list args;
+ fprintf(stderr," Diagnostic: ");
+ va_start(args, diag);
+ vfprintf(stderr, diag, args);
+ va_end(args);
+ fprintf(stderr,"\n");
+ }
+ fprintf(stderr,"usage: %s [-options] [prevcal] inoutfile\n",error_program);
+ fprintf(stderr," -v verbosity Verbose mode\n");
+ fprintf(stderr," -p Plot graphs.\n");
+ fprintf(stderr," -i Initial calibration, set targets, create .cal\n");
+ fprintf(stderr," -r Re-calibrate against previous .cal and create new .cal\n");
+ fprintf(stderr," -e Verify against previous .cal\n");
+ fprintf(stderr," -I Create imitation target from .ti3 and null calibration\n");
+ fprintf(stderr," -d Go through the motions but don't write any files\n");
+ fprintf(stderr," -s smoothing Extra curve smoothing (default 1.0)\n");
+ fprintf(stderr," -A manufacturer Set the manufacturer description string\n");
+ fprintf(stderr," -M model Set the model description string\n");
+ fprintf(stderr," -D description Set the profile Description string\n");
+ fprintf(stderr," -C copyright Set the copyright string\n");
+ fprintf(stderr," -x# percent Set initial maximum device %% target (override auto)\n");
+ fprintf(stderr," -m# percent Set initial dev target to %% of auto maximum\n");
+ fprintf(stderr," -n# deltaE Set initial white minimum deltaE target\n");
+ fprintf(stderr," -t# percent Set initial 50%% transfer curve percentage target\n");
+ fprintf(stderr," # = c, r, 0 First channel\n");
+ fprintf(stderr," m, g, 1 Second channel\n");
+ fprintf(stderr," y, b, 2 Third channel\n");
+ fprintf(stderr," k, 3 Fourth channel, etc.\n");
+ fprintf(stderr," -a Create an Adobe Photoshop .AMP file as well as a .cal\n");
+ fprintf(stderr," prevcal Base name of previous .cal file for recal or verify.\n");
+ fprintf(stderr," inoutname Base name of input .ti3 file, output .cal file\n");
+ exit(1);
+}
+
+/* - - - - - - - - - - - - - - - - - - - - - - - */
+typedef struct {
+ double loc; /* Location up the curve 0.0 - 1.0 */
+ double val[MAX_CHAN]; /* Value at that location 0.0 - 1.0 */
+} trans_point;
+
+/* Class to hold a print calibration target */
+struct _pcaltarg {
+ inkmask devmask; /* ICX ink mask of device space */
+
+ /* Note that with all of these, a value < 0.0 */
+ /* indicates no value set. */
+ int devmaxset; /* Flag - nz if the devmax is set */
+ double devmax[MAX_CHAN]; /* Device value maximum 0.0 - 1.0 */
+
+ int ademaxset; /* Flag - nz if the ademax is set */
+ double ademax[MAX_CHAN]; /* abs DE maximum for each channel */
+
+ int ademinset; /* Flag - nz if the ademin is set */
+ double ademin[MAX_CHAN]; /* abs DE minimum for each channel */
+
+ int no_tpoints; /* Number of transfer curve points */
+ trans_point *tpoints; /* Array of transfer curve points */
+
+ char err[500]; /* Error message from diagnostics */
+
+ /* Methods */
+ void (*del)(struct _pcaltarg *p);
+
+ /* Save/restore to a CGATS file */
+ int (*write)(struct _pcaltarg *p, cgats *cg, int tab); /* return nz on error */
+ int (*read)(struct _pcaltarg *p, cgats *cg, int tab); /* return nz on error */
+
+ /* Set values in the target */
+ void (*update_devmax)(struct _pcaltarg *p, int chan, double val);
+ void (*update_ademax)(struct _pcaltarg *p, int chan, double val);
+ void (*update_ademin)(struct _pcaltarg *p, int chan, double val);
+ void (*update_tcurve)(struct _pcaltarg *p, int chan, double loc, double val);
+
+ /* Reurn nz if the target has been set */
+ int (*is_set)(struct _pcaltarg *p);
+
+ /* Update settings or from one from another */
+ void (*update)(struct _pcaltarg *p, struct _pcaltarg *s);
+
+}; typedef struct _pcaltarg pcaltarg;
+
+static void pcaltarg_del(pcaltarg *p) {
+ if (p != NULL) {
+ free(p);
+ }
+}
+
+/* Write the cal target to a givent cgats table */
+static int pcaltarg_write(pcaltarg *p, cgats *cg, int tab) {
+ int i, j;
+ time_t clk = time(0);
+ struct tm *tsp = localtime(&clk);
+ char *atm = asctime(tsp); /* Ascii time */
+ char *ident = icx_inkmask2char(p->devmask, 1);
+ char *bident = icx_inkmask2char(p->devmask, 0);
+ int devchan = icx_noofinks(p->devmask);
+ 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 */
+ cg->add_table(cg, tt_other, 0); /* Add a table for Calibration TarGet values */
+ cg->add_kword(cg, tab, "DESCRIPTOR", "Argyll Calibration Target Definition File",NULL);
+ cg->add_kword(cg, tab, "ORIGINATOR", "Argyll printcal", NULL);
+ cg->add_kword(cg, tab, "CREATED",atm, NULL);
+ cg->add_kword(cg, tab, "COLOR_REP", ident, NULL);
+
+ /* Setup the table, which holds all the model parameters. */
+ /* There is always a parameter per X Y Z or spectral band */
+ cg->add_field(cg, tab, "PARAMTYPE", nqcs_t);
+ nsetel++;
+ sprintf(buf, "%s_I",bident);
+ cg->add_field(cg, tab, buf, r_t);
+ nsetel++;
+ for (i = 0; i < devchan; i++) {
+ inkmask imask = icx_index2ink(p->devmask, i);
+ sprintf(buf, "%s_%s",bident,icx_ink2char(imask));
+ cg->add_field(cg, tab, buf, r_t);
+ nsetel++;
+ }
+
+ if ((setel = (cgats_set_elem *)malloc(sizeof(cgats_set_elem) * nsetel)) == NULL) {
+ free(ident);
+ free(bident);
+ sprintf(p->err,"ctg_write: malloc of setel failed");
+ return 1;
+ }
+
+ /* Write out the values */
+ if (p->devmaxset) {
+ /* This is informational only */
+ setel[0].c = "DEVMAX_USED";
+ setel[1].d = 0.0; /* Not used */
+
+ if (p->devmask & ICX_ADDITIVE) {
+ for (i = 0; i < devchan; i++)
+ setel[2+i].d = 1.0 - p->devmax[i];
+ } else {
+ for (i = 0; i < devchan; i++)
+ setel[2+i].d = p->devmax[i];
+ }
+ cg->add_setarr(cg, tab, setel);
+ }
+ if (p->ademaxset) {
+ setel[0].c = "DELMAX_AIM";
+ setel[1].d = 0.0; /* Not used */
+
+ for (i = 0; i < devchan; i++)
+ setel[2+i].d = p->ademax[i];
+ cg->add_setarr(cg, tab, setel);
+ }
+ if (p->ademinset) {
+ setel[0].c = "DELMIN_AIM";
+ setel[1].d = 0.0; /* Not used */
+
+ for (i = 0; i < devchan; i++)
+ setel[2+i].d = p->ademin[i];
+ cg->add_setarr(cg, tab, setel);
+ }
+ for (j = 0; j < p->no_tpoints; j++) {
+ setel[0].c = "TRANS_PNT";
+ setel[1].d = p->tpoints[j].loc;
+
+ for (i = 0; i < devchan; i++)
+ setel[2+i].d = p->tpoints[j].val[i];
+ cg->add_setarr(cg, tab, setel);
+ }
+ free(setel);
+ free(ident);
+ free(bident);
+
+ return 0;
+}
+
+/* Read the cal target from a given cgats table */
+static int pcaltarg_read(pcaltarg *p, cgats *cg, int tab) {
+ char *bident;
+ int devchan;
+ int i, j, ix;
+ int ti; /* Temporary CGATs index */
+ int spi[2+MAX_CHAN]; /* CGATS indexes for each field */
+ char buf[100];
+
+ if ((ti = cg->find_kword(cg, tab, "COLOR_REP")) < 0) {
+ sprintf(p->err, "ctg_read: Can't fint COLOR_REP");
+ return 1;
+ }
+
+ if ((p->devmask = icx_char2inkmask(cg->t[tab].kdata[ti]) ) == 0) {
+ sprintf(p->err, "ctg_read: unrecognized COLOR_REP '%s'",cg->t[tab].kdata[ti]);
+ return 1;
+ }
+ devchan = icx_noofinks(p->devmask);
+ bident = icx_inkmask2char(p->devmask, 0);
+
+ /* Figure out the indexes of all the fields */
+ if ((spi[0] = cg->find_field(cg, tab, "PARAMTYPE")) < 0) {
+ sprintf(p->err, "ctg_read: Can't find field PARAMTYPE");
+ free(bident);
+ return 1;
+ }
+ sprintf(buf, "%s_I",bident);
+ if ((spi[1] = cg->find_field(cg, tab, buf)) < 0) {
+ sprintf(p->err, "ctg_read: Can't find field %s",buf);
+ free(bident);
+ return 1;
+ }
+ for (i = 0; i < devchan; i++) {
+ inkmask imask = icx_index2ink(p->devmask, i);
+ sprintf(buf, "%s_%s",bident,icx_ink2char(imask));
+ if ((spi[2+i] = cg->find_field(cg, tab, buf)) < 0) {
+ sprintf(p->err, "ctg_read: Can't find field %s",buf);
+ free(bident);
+ return 1;
+ }
+ }
+
+ /* Go through all the entries in the table, putting them in the right place */
+ for (ix = 0; ix < cg->t[tab].nsets; ix++) {
+
+ if (strcmp((char *)cg->t[tab].fdata[ix][spi[0]], "DELMAX_AIM") == 0) {
+ for (i = 0; i < devchan; i++)
+ p->ademax[i] = *((double *)cg->t[tab].fdata[ix][spi[2+i]]);
+ p->ademaxset = 1;
+
+ } else if (strcmp((char *)cg->t[tab].fdata[ix][spi[0]], "DELMIN_AIM") == 0) {
+ for (i = 0; i < devchan; i++)
+ p->ademin[i] = *((double *)cg->t[tab].fdata[ix][spi[2+i]]);
+ p->ademinset = 1;
+
+ } else if (strcmp((char *)cg->t[tab].fdata[ix][spi[0]], "TRANS_PNT") == 0) {
+ if ((p->tpoints = (trans_point *)realloc(p->tpoints, sizeof(trans_point)
+ * (p->no_tpoints+1))) == NULL)
+ error("Realloc of tpoints");
+ p->tpoints[p->no_tpoints].loc = *((double *)cg->t[tab].fdata[ix][spi[1]]);
+ for (i = 0; i < devchan; i++)
+ p->tpoints[p->no_tpoints].val[i] = *((double *)cg->t[tab].fdata[ix][spi[2+i]]);
+ p->no_tpoints++;
+ }
+ }
+ free(bident);
+
+ return 0;
+}
+
+/* Update an individual setting. Use chan < 0 to set all to default */
+void pcaltarg_update_devmax(struct _pcaltarg *p, int chan, double val) {
+ int i;
+ if (p->devmaxset == 0) {
+ for (i = 0; i < MAX_CHAN; i++)
+ p->devmax[i] = -1.0;
+ p->devmaxset = 1;
+ }
+ if (chan >= 0)
+ p->devmax[chan] = val;
+}
+void pcaltarg_update_ademax(struct _pcaltarg *p, int chan, double val) {
+ int i;
+ if (p->ademaxset == 0) {
+ for (i = 0; i < MAX_CHAN; i++)
+ p->ademax[i] = -1.0;
+ p->ademaxset = 1;
+ }
+ if (chan >= 0)
+ p->ademax[chan] = val;
+}
+void pcaltarg_update_ademin(struct _pcaltarg *p, int chan, double val) {
+ int i;
+ if (p->ademinset == 0) {
+ for (i = 0; i < MAX_CHAN; i++)
+ p->ademin[i] = -1.0;
+ p->ademinset = 1;
+ }
+ if (chan >= 0)
+ p->ademin[chan] = val;
+}
+void pcaltarg_update_tcurve(struct _pcaltarg *p, int chan, double loc, double val) {
+ int i, j;
+
+ /* See if a transfer curve point already exists */
+ for (j = 0; j < p->no_tpoints; j++) {
+ if (p->tpoints[j].loc == loc)
+ break;
+ }
+ /* If not, allocate a new one */
+ if (j >= p->no_tpoints) {
+ p->no_tpoints++;
+ if ((p->tpoints = (trans_point *)realloc(p->tpoints, sizeof(trans_point) * p->no_tpoints)) == NULL)
+ error("Realloc of tpoints");
+ p->tpoints[j].loc = loc;
+ for (i = 0; i < MAX_CHAN; i++)
+ p->tpoints[j].val[i] = -1.0;
+ }
+ p->tpoints[j].val[chan] = val;
+
+ if (p->no_tpoints > 0) {
+ /* Sort the transfer points into loc order */
+#define HEAP_COMPARE(A,B) ((A).loc < (B).loc)
+ HEAPSORT(trans_point, p->tpoints, p->no_tpoints);
+#undef HEAP_COMPARE
+ }
+}
+
+/* Reurn nz if the target has been set */
+static int pcaltarg_is_set(pcaltarg *p) {
+ if (p->devmaxset != 0
+ || p->ademaxset != 0
+ || p->ademinset != 0
+ || p->no_tpoints > 0)
+ return 1;
+ return 0;
+}
+
+/* Update one from another */
+static void pcaltarg_update(pcaltarg *p, pcaltarg *s) {
+ int i, j, k;
+
+ if (s->devmaxset) {
+ if (p->devmaxset == 0) {
+ for (i = 0; i < MAX_CHAN; i++)
+ p->devmax[i] = -1.0;
+ p->devmaxset = 1;
+ }
+ for (i = 0; i < MAX_CHAN; i++) {
+ if (s->devmax[i] >= 0.0)
+ p->devmax[i] = s->devmax[i];
+ }
+ }
+
+ if (s->ademaxset) {
+ if (p->ademaxset == 0) {
+ for (i = 0; i < MAX_CHAN; i++)
+ p->ademax[i] = -1.0;
+ p->ademaxset = 1;
+ }
+ for (i = 0; i < MAX_CHAN; i++) {
+ if (s->ademax[i] >= 0.0)
+ p->ademax[i] = s->ademax[i];
+ }
+ }
+
+ if (s->ademinset) {
+ if (p->ademinset == 0) {
+ for (i = 0; i < MAX_CHAN; i++)
+ p->ademin[i] = -1.0;
+ p->ademinset = 1;
+ }
+ for (i = 0; i < MAX_CHAN; i++) {
+ if (s->ademin[i] >= 0.0)
+ p->ademin[i] = s->ademin[i];
+ }
+ }
+
+ /* For each source transfer curve point */
+ for (k = 0; k < s->no_tpoints; k++) {
+
+ /* See if a transfer curve point already exists */
+ for (j = 0; j < p->no_tpoints; j++) {
+ if (p->tpoints[j].loc == s->tpoints[k].loc)
+ break;
+ }
+ /* If not, allocate a new one */
+ if (j >= p->no_tpoints) {
+ p->no_tpoints++;
+ if ((p->tpoints = (trans_point *)realloc(p->tpoints, sizeof(trans_point) * p->no_tpoints)) == NULL)
+ error("Realloc of tpoints");
+ p->tpoints[j].loc = s->tpoints[k].loc;
+ for (i = 0; i < MAX_CHAN; i++)
+ p->tpoints[j].val[i] = -1.0;
+ }
+ for (i = 0; i < MAX_CHAN; i++) {
+ if (s->tpoints[k].val[i] >= 0.0)
+ p->tpoints[j].val[i] = s->tpoints[k].val[i];
+ }
+ }
+ if (s->no_tpoints > 0) {
+ /* Sort the transfer points into loc order */
+#define HEAP_COMPARE(A,B) ((A).loc < (B).loc)
+ HEAPSORT(trans_point, p->tpoints, p->no_tpoints);
+#undef HEAP_COMPARE
+ }
+}
+
+/* Create a new, empty pcaltarget */
+/* Return NULL on error */
+pcaltarg *new_pcaltarg() {
+ pcaltarg *p;
+
+ if ((p = (pcaltarg *)calloc(1, sizeof(pcaltarg))) == NULL) {
+ return NULL;
+ }
+
+ /* Set method pointers */
+ p->del = pcaltarg_del;
+ p->write = pcaltarg_write;
+ p->read = pcaltarg_read;
+ p->update_devmax = pcaltarg_update_devmax;
+ p->update_ademax = pcaltarg_update_ademax;
+ p->update_ademin = pcaltarg_update_ademin;
+ p->update_tcurve = pcaltarg_update_tcurve;
+ p->is_set = pcaltarg_is_set;
+ p->update = pcaltarg_update;
+
+ return p;
+}
+
+/* - - - - - - - - - - - - - - - - - - - - - - - */
+
+/* A wedge sample value */
+typedef struct {
+ double inv; /* Input value (cal table) */
+ double dev; /* Device value */
+ double XYZ[3]; /* XYZ value */
+ double Lab[3]; /* Lab value */
+ double del; /* Absolute delta (to white) */
+} wval;
+
+#define MAX_INVSOLN 10 /* Rspl maximum reverse solutions */
+
+/* rspl setting functions */
+static void rsplset1(void *cbntx, double *out, double *in) {
+ co *dpoints = (co *)cbntx;
+ int ix;
+
+ ix = *((int*)&in[-0-1]); /* Get grid index being looked up */
+ out[0] = dpoints[ix].v[0];
+}
+
+/* Do an inverse lookup of an rspl. Return -1.0 on error. */
+/* dir is value to favour if there are multiple solutions. */
+static double rspl_ilookup(rspl *r, double dir, double in) {
+ int nsoln; /* Number of solutions found */
+ co pp[MAX_INVSOLN]; /* Room for all the solutions found */
+ int k; /* Chosen solution */
+
+ pp[0].v[0] = in;
+
+ nsoln = r->rev_interp (
+ r, /* this */
+ RSPL_NEARCLIP, /* Clip to nearest (faster than vector) */
+ MAX_INVSOLN, /* Maximum number of solutions allowed for */
+ NULL, /* No auxiliary input targets */
+ NULL, /* Clip vector direction and length */
+ pp); /* Input and output values */
+
+ nsoln &= RSPL_NOSOLNS; /* Get number of solutions */
+
+ if (nsoln == 1) { /* Exactly one solution */
+ k = 0;
+ } else if (nsoln == 0) { /* Zero solutions. This is unexpected. */
+ return -1.0;
+ } else { /* Multiple solutions */
+ double bdist = 1e300;
+ int bsoln = 0;
+// warning("Multiple solutions for curve %d for DE %f",j,pp[0].v[0]);
+ for (k = 0; k < nsoln; k++) {
+ double tt;
+ tt = pp[k].p[0] - dir;
+ tt *= tt;
+ if (tt < bdist) { /* Better solution */
+ bdist = tt;
+ bsoln = k;
+ }
+ }
+ k = bsoln;
+ }
+ return pp[k].p[0];
+}
+
+int main(int argc, char *argv[]) {
+ int fa,nfa,mfa; /* current argument we're looking at */
+ int verb = 0;
+ int doplot = 0;
+ int initial = 0; /* Do initial creation of cal target and calibration */
+ int recal = 0; /* Do recalibrate/use cal target. */
+ int verify = 0; /* Do verification */
+ int imitate = 0; /* Do target directly from input */
+ int dowrite = 1; /* Write to files */
+ int doamp = 0; /* Write Adobe Photoshop .AMP file */
+ profxinf xpi; /* Extra profile/calibration information */
+ pcaltarg *upct = NULL; /* User settings of print calibration target */
+ pcaltarg *pct = NULL; /* Settings of print calibration target */
+ double smooth = RSPLSMOOTH; /* RSPL Smoothness factor */
+ double xsmooth = 1.0; /* Smoothing multiplier */
+ double ver_maxde = 2.0; /* Verify maximum Delta E (1.0 for smooth == 1.0) */
+ int spec = 0; /* Use spectral data flag */
+ icxIllumeType illum = icxIT_D50; /* Spectral defaults */
+ xspect cust_illum; /* Custom illumination spectrum */
+ icxObserverType observ = icxOT_CIE_1931_2; /* The classic observer */
+
+ char baname[MAXNAMEL+1] = ""; /* Input & Output base name */
+ char inname[MAXNAMEL+1] = ""; /* new .ti3 input file name */
+ char calname[MAXNAMEL+1] = ""; /* previous .cal input file name */
+ char outname[MAXNAMEL+1] = ""; /* new .cal output file name */
+ char ampname[MAXNAMEL+1] = ""; /* new .amp output file name */
+ double maxscale[MAX_CHAN]; /* Scale auto device maximum to % */
+ cgats *icg = NULL; /* .ti3 input cgats structure */
+ int ti; /* Temporary CGATs index */
+ inkmask devmask; /* ICX ink mask of device space */
+ int devchan; /* Number of chanels in device space */
+ int isLab = 0; /* Flag indicating whether PCS is XYZ or Lab */
+ int n_pvals[MAX_CHAN]; /* Number of measurement values */
+ wval *pvals[MAX_CHAN]; /* Patch measurement values */
+ wval white; /* Average white value */
+ int n_white = 0; /* Number of values to average */
+ icmXYZNumber wht; /* White value */
+ rspl *raw[MAX_CHAN]; /* Raw Lab values fitted to rspl */
+ rspl *ade[MAX_CHAN]; /* Absolute delta E */
+ rspl *rde[MAX_CHAN]; /* Relative delta E */
+ rspl *pcade[MAX_CHAN]; /* Previous calibrated absolute delta E */
+ double mxade[MAX_CHAN]; /* Maximum ade value */
+ double idpow[MAX_CHAN] = { -1.0 }; /* Ideal power-like of targen values */
+ int n_cvals; /* Number of calibration curve values */
+ wval *cvals[MAX_CHAN]; /* Calibration curve tables */
+ rspl *tcurves[MAX_CHAN]; /* Tweak target curves */
+
+ int i, j;
+
+ /* Init pointers to NULL */
+ for (j = 0; j < MAX_CHAN; j++) {
+ maxscale[j] = -1.0;
+ pvals[j] = NULL;
+ raw[j] = NULL;
+ ade[j] = NULL;
+ rde[j] = NULL;
+ pcade[j] = NULL;
+ cvals[j] = NULL;
+ tcurves[j] = NULL;
+ }
+
+ error_program = argv[0];
+ memset((void *)&xpi, 0, sizeof(profxinf)); /* Init extra profile info to defaults */
+ if ((upct = new_pcaltarg()) == NULL || (pct = new_pcaltarg()) == NULL)
+ error("new_caltarg failed");
+
+ if (argc < 3)
+ usage("Too few arguments, got %d expect at least %d",argc-1,2);
+
+#ifdef NEVER
+ {
+ double src, dst, pp;
+
+ src = 0.5;
+ dst = 0.25;
+ pp = icx_powlike_needed(src, dst);
+ printf("%f -> %f needs %f, check %f\n",src,dst,pp,icx_powlike(src,pp));
+
+ src = 0.25;
+ dst = 0.5;
+ pp = icx_powlike_needed(src, dst);
+ printf("%f -> %f needs %f, check %f\n",src,dst,pp,icx_powlike(src,pp));
+
+ src = 0.5;
+ dst = 0.707106;
+ pp = icx_powlike_needed(src, dst);
+ printf("%f -> %f needs %f, check %f\n",src,dst,pp,icx_powlike(src,pp));
+
+ src = 0.5;
+ dst = 0.5;
+ pp = icx_powlike_needed(src, dst);
+ printf("%f -> %f needs %f, check %f\n",src,dst,pp,icx_powlike(src,pp));
+ }
+#endif // NEVER
+
+ /* Process the arguments */
+ mfa = 2; /* Minimum final arguments */
+ for(fa = 1;fa < argc;fa++) {
+
+ nfa = fa; /* skip to nfa if next argument is used */
+ if (argv[fa][0] == '-') { /* Look for any flags */
+ char *na = NULL; /* next argument after flag, null if none */
+
+ if (argv[fa][2] != '\000')
+ na = &argv[fa][2]; /* next is directly after flag */
+ else {
+ if ((fa+1+mfa) < argc) {
+ if (argv[fa+1][0] != '-') {
+ nfa = fa + 1;
+ na = argv[nfa]; /* next is seperate non-flag argument */
+ }
+ }
+ }
+
+ if (argv[fa][1] == '?')
+ usage("Usage requested");
+
+ else if (argv[fa][1] == 'v') {
+ if (na != NULL) {
+ fa = nfa;
+ verb = atoi(na);
+ } else
+ verb = 1;
+ }
+
+ else if (argv[fa][1] == 'p') {
+ if (na != NULL) {
+ fa = nfa;
+ doplot = atoi(na);
+ } else
+ doplot = 1; /* Plot various graphs */
+ }
+
+ else if (argv[fa][1] == 'i') {
+ initial = 1; /* Initial calibration */
+ recal = 0;
+ verify = 0;
+ imitate = 0;
+ }
+
+ else if (argv[fa][1] == 'r') {
+ initial = 0;
+ recal = 1; /* Recalibrate */
+ verify = 0;
+ imitate = 0;
+ }
+
+ else if (argv[fa][1] == 'e') {
+ initial = 0;
+ recal = 0;
+ verify = 1; /* Verify */
+ imitate = 0;
+ }
+
+ else if (argv[fa][1] == 'I') {
+ initial = 0;
+ recal = 0;
+ verify = 0;
+ imitate = 1; /* Imitation target */
+ }
+
+ else if (argv[fa][1] == 'd')
+ dowrite = 0; /* Don't write to files */
+
+ else if (argv[fa][1] == 'a')
+ doamp = 1; /* write AMP file */
+
+ /* Smoothing modfider */
+ else if (argv[fa][1] == 's') {
+ fa = nfa;
+ if (na == NULL) usage("Expect argument to smoothing flag -s");
+ xsmooth = atof(na);
+ }
+
+ /* Manufacturer description string */
+ else if (argv[fa][1] == 'A') {
+ fa = nfa;
+ if (na == NULL) usage("Expect argument to manufacturer description flag -A");
+ xpi.deviceMfgDesc = na;
+ }
+
+ /* Model description string */
+ else if (argv[fa][1] == 'M') {
+ fa = nfa;
+ if (na == NULL) usage("Expect argument to model description flag -M");
+ xpi.modelDesc = na;
+ }
+
+ /* Profile Description */
+ else if (argv[fa][1] == 'D') {
+ fa = nfa;
+ if (na == NULL) usage("Expect argument to profile description flag -D");
+ xpi.profDesc = na;
+ }
+
+ /* Copyright string */
+ else if (argv[fa][1] == 'C') {
+ fa = nfa;
+ if (na == NULL) usage("Expect argument to copyright flag -C");
+ xpi.copyright = na;
+ }
+
+ /* Per channel target modifiers */
+ else if (argv[fa][1] == 'x'
+ || argv[fa][1] == 'm'
+ || argv[fa][1] == 'n'
+ || argv[fa][1] == 't') {
+ char fch = argv[fa][1];
+ int chan = -1;
+ double val = -1.0;
+ fa = nfa;
+ if (na == NULL)
+ usage("Expect channel flag after flag -%c",argv[fa][1]);
+ switch (na[0]) {
+ case 'c': case 'r': case '0':
+ chan = 0;
+ break;
+ case 'm': case 'g': case '1':
+ chan = 1;
+ break;
+ case 'y': case 'b': case '2':
+ chan = 2;
+ break;
+ case 'k': case '3':
+ chan = 3;
+ break;
+ case '4':
+ chan = 4;
+ break;
+ case '5':
+ chan = 5;
+ break;
+ case '6':
+ chan = 6;
+ break;
+ case '7':
+ chan = 7;
+ break;
+ case '8':
+ chan = 8;
+ break;
+ case '9':
+ chan = 9;
+ break;
+ case 'A':
+ chan = 10;
+ break;
+ case 'B':
+ chan = 11;
+ break;
+ case 'C':
+ chan = 12;
+ break;
+ case 'D':
+ chan = 13;
+ break;
+ case 'E':
+ chan = 14;
+ break;
+ case 'F':
+ chan = 15;
+ break;
+ default:
+ usage("Unknown channel flag '%s' after flag -%c",argv[fa][2],argv[fa][1]);
+ }
+ ++fa;
+ if (fa >= argc || argv[fa][0] == '-') usage("Expect argument after flag -%c%c",fch,na[0]);
+ val = atof(argv[fa]);
+
+ if (fch == 'x') {
+ if (val < 0.0 || val > 100.0)
+ usage("Argument to -%c%c %f from '%s' is out of range",fch,na[0],val,argv[fa]);
+ val /= 100.0;
+ upct->update_devmax(upct, chan, val);
+
+ } else if (fch == 'm') {
+ if (val < 0.0 || val > 100.0)
+ usage("Argument to -%c%c %f from '%s' is out of range",fch,na[0],val,argv[fa]);
+ val /= 100.0;
+ maxscale[chan] = val;
+
+ } else if (fch == 'n') {
+ upct->update_ademin(upct, chan, val);
+
+ } else if (fch == 't') {
+ if (val < 0.0 || val > 100.0)
+ usage("Argument to -%c%c %f from '%s' is out of range",fch,na[0],val,argv[fa]);
+ val /= 100.0;
+ upct->update_tcurve(upct, chan, 0.5, val);
+ }
+ }
+ else
+ usage("Unknown flag '%c'",argv[fa][1]);
+ } else
+ break;
+ }
+
+ smooth *= xsmooth;
+
+ if (!( (initial && !recal && !verify && !imitate)
+ || (!initial && recal && !verify && !imitate)
+ || (!initial && !recal && verify && !imitate)
+ || (!initial && !recal && !verify && imitate)))
+ error("One of -i, -r -e or -I must be set");
+
+ /* Get the file name arguments */
+ if (verify || recal) {
+ if (fa >= argc || argv[fa][0] == '-') usage("Missing prevoius .cal basename");
+ strncpy(calname,argv[fa++],MAXNAMEL-4); calname[MAXNAMEL-4] = '\000';
+ strcat(calname,".cal");
+ }
+ if (fa >= argc || argv[fa][0] == '-') usage("Missing .ti3 and new .cal basename");
+ strncpy(baname,argv[fa++],MAXNAMEL-4); baname[MAXNAMEL-4] = '\000';
+ strcpy(inname,baname); /* new .ti3 file */
+ strcat(inname,".ti3");
+ strcpy(outname,baname); /* New .cal file */
+ strcat(outname,".cal");
+ strcpy(ampname,baname); /* New .amp file */
+ strcat(ampname,".amp");
+
+ if (fa < argc) usage("Too many arguments ('%s')",argv[fa]);
+
+ /* Open and look at the .ti3 profile patches file */
+ icg = new_cgats(); /* Create a CGATS structure */
+ icg->add_other(icg, "CTI3"); /* our special input type is Calibration Target Information 3 */
+
+ if (icg->read_name(icg, inname))
+ error("CGATS file read error : %s",icg->err);
+
+ if (icg->ntables == 0 || icg->t[0].tt != tt_other || icg->t[0].oi != 0)
+ error ("Input file isn't a CTI3 format file");
+ if (icg->ntables < 1)
+ error ("Input file doesn't contain at least one table");
+
+ /* See if CIE is actually available - some sources of .TI3 don't provide it */
+ if (!spec
+ && icg->find_field(icg, 0, "LAB_L") < 0
+ && icg->find_field(icg, 0, "XYZ_X") < 0) {
+
+ if (icg->find_kword(icg, 0, "SPECTRAL_BANDS") < 0)
+ error ("Neither CIE nor spectral data found in file '%s'",inname);
+
+ /* Switch to using spectral information */
+ if (verb)
+ printf("No CIE data found, switching to spectral with standard observer & D50\n");
+ spec = 1;
+ illum = icxIT_D50;
+ observ = icxOT_CIE_1931_2;
+ }
+
+ /* If we requested spectral, check that it is available */
+ if (spec) {
+ if (icg->find_kword(icg, 0, "SPECTRAL_BANDS") < 0)
+ error ("Requested spectral interpretation when data not available");
+ }
+
+ /* Get colorspace information from input CGATS file */
+ {
+ char *buf;
+ char *inc, *outc;
+
+ if ((ti = icg->find_kword(icg, 0, "COLOR_REP")) < 0)
+ error("Input file doesn't contain keyword COLOR_REPS");
+
+ if ((buf = strdup(icg->t[0].kdata[ti])) == NULL)
+ error("Malloc failed - color rep");
+
+ /* Split COLOR_REP into device and PCS space */
+ inc = buf;
+ if ((outc = strchr(buf, '_')) == NULL)
+ error("COLOR_REP '%s' invalid", icg->t[0].kdata[ti]);
+ *outc++ = '\000';
+
+ if (strcmp(outc, "XYZ") == 0)
+ isLab = 0;
+ else if (strcmp(outc, "LAB") == 0)
+ isLab = 1;
+ else
+ error("COLOR_REP '%s' invalid (Neither XYZ nor LAB)", icg->t[0].kdata[ti]);
+
+ devmask = icx_char2inkmask(inc);
+ devchan = icx_noofinks(devmask);
+
+ if (devchan == 0)
+ error("COLOR_REP '%s' invalid (No matching devmask)", icg->t[0].kdata[ti]);
+
+ if ((devmask & ICX_ADDITIVE) && !(devmask & ICX_INVERTED))
+ warning("COLOR_REP '%s' is probably not suitable for print calibration!", icg->t[0].kdata[ti]);
+
+ free(buf);
+ }
+
+ if (verify || recal || imitate) {
+ if (upct->is_set(upct)) {
+ warning("Command line calibration target paramers ignored on re-calibrate, verify and imitate!");
+ }
+ }
+
+ /* For recalibrate or verify, load the previous calibration file */
+ if (verify || recal) {
+ cgats *tcg; /* Previous .cal file */
+
+ tcg = new_cgats(); /* Create a CGATS structure */
+ tcg->add_other(tcg, "CAL"); /* our special input type is Calibration Target */
+
+ if (tcg->read_name(tcg, calname))
+ error("No cal target '%s' found for re-calibrate (%s)\n",calname,tcg->err);
+
+ /* Check that this is an output cal file */
+ if ((ti = tcg->find_kword(tcg, 0, "DEVICE_CLASS")) < 0)
+ error ("Calibration file '%s'doesn't contain keyword DEVICE_CLASS",calname);
+ if (strcmp(tcg->t[0].kdata[ti],"OUTPUT") != 0)
+ error ("Calibration file '%s' doesn't has DEVICE_CLASS that is not OUTPUT",calname);
+
+ if (pct->read(pct, tcg, 1) != 0)
+ error("Reading cal target '%s' failed",calname);
+
+ if (pct->devmask != devmask)
+ error("Target '%s' colorspace '%s' doesn't match '%s' colorspace '%s'",
+ calname,icx_inkmask2char(pct->devmask, 1),inname,icx_inkmask2char(devmask, 1));
+
+ /* Load the previous expected absolute DE response */
+ /* It will be in the third table with other type "CAL" */
+ if (tcg->ntables >= 3 && tcg->t[2].tt == tt_other && tcg->t[0].oi == 0) {
+ int ti;
+ char *bident;
+ int spi[1+MAX_CHAN]; /* CGATS indexes for each field */
+ char buf[100];
+
+ bident = icx_inkmask2char(pct->devmask, 0);
+
+ if (tcg->t[2].nsets <= 0)
+ error ("No Calibration Expected DE Response in '%s'",calname);
+
+ /* Figure out the indexes of all the fields */
+ sprintf(buf, "%s_I_DE",bident);
+ if ((spi[0] = tcg->find_field(tcg, 2, buf)) < 0)
+ error("Can't find field %s in '%s'",buf,calname);
+
+ for (i = 0; i < devchan; i++) {
+ inkmask imask = icx_index2ink(pct->devmask, i);
+ sprintf(buf, "%s_%s_DE",bident,icx_ink2char(imask));
+ if ((spi[1+i] = tcg->find_field(tcg, 2, buf)) < 0)
+ error("Can't find field %s in '%s'",buf,calname);
+ }
+
+ /* Read in each channels values and put them in a rspl */
+ for (j = 0; j < devchan; j++) {
+ datai low,high;
+ int gres[MXDI];
+ co *dpoints;
+
+ low[0] = 0.0;
+ high[0] = 1.0;
+ gres[0] = tcg->t[2].nsets;
+
+ if ((pcade[j] = new_rspl(RSPL_NOFLAGS,1, 1)) == NULL)
+ error("new_rspl() failed");
+
+ if ((dpoints = malloc(sizeof(co) * gres[0])) == NULL)
+ error("malloc dpoints[%d] failed",gres[0]);
+
+ /* Copy the points to our array */
+ if (devmask & ICX_ADDITIVE) {
+ for (i = 0; i < gres[0]; i++) {
+ dpoints[i].p[0] = 1.0 - i/(double)(gres[0]-1);
+ dpoints[i].v[0] = *((double *)tcg->t[2].fdata[gres[0]-1-i][spi[1+j]]);
+ }
+ } else {
+ for (i = 0; i < gres[0]; i++) {
+ dpoints[i].p[0] = i/(double)(gres[0]-1);
+ dpoints[i].v[0] = *((double *)tcg->t[2].fdata[i][spi[1+j]]);
+ }
+ }
+
+ pcade[j]->set_rspl(pcade[j],
+ 0,
+ (void *)dpoints, /* Read points */
+ rsplset1, /* Setting function */
+ low, high, gres, /* Low, high, resolution of grid */
+ NULL, NULL /* Default data scale */
+ );
+ free(dpoints);
+ }
+ free(bident);
+ }
+ tcg->del(tcg);
+
+ } else { /* Must be an initial or Imitation calibration */
+
+ pct->devmask = devmask;
+
+ /* Set the cal target from any user supplied parameters */
+ pct->update(pct, upct);
+
+ /* No previous absolute de reference */
+ for (j = 0; j < devchan; j++)
+ pcade[j] = NULL;
+ }
+
+ /* Common processing: */
+
+ /* Read in the patch data */
+ {
+ char buf[100];
+ char *pcsfname[2][3] = { { "XYZ_X", "XYZ_Y", "XYZ_Z" },
+ { "LAB_L", "LAB_A", "LAB_B" } };
+ int dvi[MAX_CHAN]; /* CGATS indexes for each device field */
+ int pcsix[3]; /* XYZ/Lab chanel indexes */
+ xsp2cie *sp2cie = NULL; /* Spectral conversion object */
+ xspect sp;
+ int spi[XSPECT_MAX_BANDS]; /* CGATS indexes for each wavelength */
+ char *bident = icx_inkmask2char(devmask, 0);
+
+ /* Figure out the indexes of all the device fields */
+ for (j = 0; j < devchan; j++) {
+ inkmask imask = icx_index2ink(devmask, j);
+ sprintf(buf, "%s_%s",bident,icx_ink2char(imask));
+ if ((dvi[j] = icg->find_field(icg, 0, buf)) < 0)
+ error("Can't find field %s in '%s'",buf,inname);
+#ifdef DEBUG
+ printf("devn chan %d field %s = %d\n",j,buf,dvi[j]);
+#endif
+ }
+ free(bident);
+
+ if (spec) {
+ int ii;
+ char buf[100];
+
+ if ((ii = icg->find_kword(icg, 0, "SPECTRAL_BANDS")) < 0)
+ error ("Input file doesn't contain keyword SPECTRAL_BANDS");
+ sp.spec_n = atoi(icg->t[0].kdata[ii]);
+ if ((ii = icg->find_kword(icg, 0, "SPECTRAL_START_NM")) < 0)
+ error ("Input file doesn't contain keyword SPECTRAL_START_NM");
+ sp.spec_wl_short = atof(icg->t[0].kdata[ii]);
+ if ((ii = icg->find_kword(icg, 0, "SPECTRAL_END_NM")) < 0)
+ error ("Input file doesn't contain keyword SPECTRAL_END_NM");
+ sp.spec_wl_long = atof(icg->t[0].kdata[ii]);
+ sp.norm = 100.0;
+
+ /* Find the fields for spectral values */
+ for (j = 0; j < sp.spec_n; j++) {
+ int nm;
+
+ /* Compute nearest integer wavelength */
+ nm = (int)(sp.spec_wl_short + ((double)j/(sp.spec_n-1.0))
+ * (sp.spec_wl_long - sp.spec_wl_short) + 0.5);
+
+ sprintf(buf,"SPEC_%03d",nm);
+
+ if ((spi[j] = icg->find_field(icg, 0, buf)) < 0)
+ error("Input file doesn't contain field %s",buf);
+ }
+
+ /* Create a spectral conversion object to XYZ */
+ if ((sp2cie = new_xsp2cie(illum, &cust_illum, observ, NULL, icSigXYZData, icxClamp)) == NULL)
+ error("Creation of spectral conversion object failed");
+
+ /* To add FWA comp. would have to locate/create spectral white here, */
+ /* then set the FWA comp. on. */
+ /* See profout.c */
+
+ } else {
+ /* Figure out the indexes of the PCS fields */
+ for (j = 0; j < 3; j++) {
+ if ((i = icg->find_field(icg, 0, pcsfname[isLab][j])) >= 0) {
+ if (icg->t[0].ftype[i] != r_t)
+ error ("Field %s is wrong type",pcsfname[isLab][j]);
+ pcsix[j] = i;
+#ifdef DEBUG
+ printf("PCS chan %d field %s = %d\n",j,pcsfname[isLab][j],pcsix[j]);
+#endif
+ } else {
+ error ("Failed to find field %s",pcsfname[isLab][j]);
+ }
+ }
+ }
+
+ n_cvals = 0;
+ for (j = 0; j < devchan; j++) {
+ pvals[j] = NULL;
+ n_pvals[j] = 0;
+ cvals[j] = NULL;
+ }
+
+ /* Read all the test patches in */
+ for (i = 0; i < icg->t[0].nsets; i++) {
+ double maxv = -1.0;
+ int maxch;
+
+#ifdef DEBUG
+ printf("Reading patch %d\n",i);
+#endif
+
+ /* Locate the maximum device value of any channel */
+ for (j = 0; j < devchan; j++) {
+ double val = *((double *)icg->t[0].fdata[i][dvi[j]]) / 100.0;
+ if (devmask & ICX_ADDITIVE)
+ val = 1.0 - val;
+ if (val > maxv) {
+ maxv = val;
+ maxch = j;
+ }
+ }
+#ifdef DEBUG
+ printf("max %f at chan %d\n",maxv,maxch);
+#endif
+ /* Treat white specially, and take it out of the list */
+ if (maxv < 1e-6) {
+ double wxyz[3];
+ if (n_white == 0) {
+ white.dev = 0.0;
+ for (j = 0; j < 3; j++)
+ white.XYZ[j] = 0.0;
+ }
+ if (spec) {
+ /* Read and convert the spectral value */
+ for (j = 0; j < sp.spec_n; j++)
+ sp.spec[j] = *((double *)icg->t[0].fdata[i][spi[j]]);
+
+ sp2cie->convert(sp2cie, wxyz, &sp);
+
+ } else {
+ /* Read the CIE value */
+ for (j = 0; j < 3; j++)
+ wxyz[j] = *((double *)icg->t[0].fdata[i][pcsix[j]]);
+
+ /* And convert to XYZ 0..1 */
+ if (isLab) {
+ icmLab2XYZ(&icmD50, wxyz, wxyz);
+ } else {
+ for (j = 0; j < 3; j++)
+ wxyz[j] /= 100.0;
+ }
+ }
+
+ white.dev += maxv;
+ for (j = 0; j < 3; j++)
+ white.XYZ[j] += wxyz[j];
+ n_white++;
+#ifdef DEBUG
+ printf(" white: dev %f,XYZ %f %f %f\n",
+ white.dev,white.XYZ[0]/n_white, white.XYZ[1]/n_white, white.XYZ[2]/n_white);
+#endif
+ } else {
+ wval *vp;
+ /* Check that all the non-max value channels are zero */
+ for (j = 0; j < devchan; j++) {
+ double val = *((double *)icg->t[0].fdata[i][dvi[j]]) / 100.0;
+ if (devmask & ICX_ADDITIVE)
+ val = 1.0 - val;
+ if (j == maxch)
+ continue;
+ if (val > 0.001)
+ break;
+ }
+ if (j < devchan) {
+#ifdef DEBUG
+ printf("Skipping patch\n");
+#endif
+ continue; /* Ignore this patch */
+ }
+
+ if ((pvals[maxch] = (wval *)realloc(pvals[maxch],
+ sizeof(wval) * (n_pvals[maxch]+1))) == NULL)
+ error("Realloc of pvals failed");
+ vp = &pvals[maxch][n_pvals[maxch]];
+ vp->dev = maxv;
+
+ if (spec) {
+ /* Read and convert the spectral value */
+ for (j = 0; j < sp.spec_n; j++)
+ sp.spec[j] = *((double *)icg->t[0].fdata[i][spi[j]]);
+
+ sp2cie->convert(sp2cie, vp->XYZ, &sp);
+
+ } else {
+ /* Read the CIE value */
+ for (j = 0; j < 3; j++)
+ vp->XYZ[j] = *((double *)icg->t[0].fdata[i][pcsix[j]]);
+
+ /* And convert to XYZ 0..1 */
+ if (isLab) {
+ icmLab2XYZ(&icmD50, vp->XYZ, vp->XYZ);
+ } else {
+ for (j = 0; j < 3; j++)
+ vp->XYZ[j] /= 100.0;
+ }
+ }
+ /* Temporary D50 Lab */
+ icmXYZ2Lab(&icmD50, vp->Lab, vp->XYZ);
+#ifdef DEBUG
+ printf(" patch %d: dev %f,XYZ %f %f %f, D50 Lab %f %f %f\n",
+ n_pvals[maxch], vp->dev,vp->XYZ[0], vp->XYZ[1], vp->XYZ[2],
+ vp->Lab[0], vp->Lab[1], vp->Lab[2]);
+#endif
+ n_pvals[maxch]++;
+ }
+ }
+
+ /* Average the white */
+ if (n_white == 0)
+ error("Can't find even one white patch in '%s'",inname);
+#ifdef DEBUG
+ printf("% white patches\n",n_white);
+#endif
+ for (j = 0; j < 3; j++) {
+ white.dev /= (double)n_white;
+ white.XYZ[j] /= (double)n_white;
+ }
+ icmAry2XYZ(wht, white.XYZ);
+ icmXYZ2Lab(&icmD50, white.Lab, white.XYZ);
+
+ /* Convert the Lab white reference to absolute */
+ wht.X /= wht.Y;
+ wht.Z /= wht.Y;
+ wht.Y /= wht.Y;
+
+ if (verb) {
+ icmXYZ2Lab(&icmD50, white.Lab, white.XYZ);
+ printf("Average white = XYZ %f %f %f, D50 Lab %f %f %f\n",
+ white.XYZ[0], white.XYZ[1], white.XYZ[2], white.Lab[0], white.Lab[1], white.Lab[2]);
+ }
+
+ for (j = 0; j < devchan; j++) {
+ wval *wp;
+
+ /* Add averaged white back into each channel */
+ if ((pvals[j] = (wval *)realloc(pvals[j],
+ sizeof(wval) * (n_pvals[j]+1))) == NULL)
+ error("Realloc (%d) of pvals failed",n_pvals[j]+1);
+ wp = &pvals[j][n_pvals[j]];
+ wp->dev = white.dev;
+ wp->XYZ[0] = white.XYZ[0];
+ wp->XYZ[1] = white.XYZ[1];
+ wp->XYZ[2] = white.XYZ[2];
+ n_pvals[j]++;
+
+ /* Convert all the XYZ values to Lab paper relative */
+ for (i = 0; i < n_pvals[j]; i++) {
+ wp = &pvals[j][i];
+ icmXYZ2Lab(&wht, wp->Lab, wp->XYZ);
+ }
+
+ /* Sort the channel acording to device value */
+ /* For a consistent result for identical device values, */
+ /* secondary sort by inverse CIE value */
+//#define HEAP_COMPARE(A,B) ((A).dev < (B).dev)
+#define HEAP_COMPARE(A,B) ((A).dev != (B).dev ? ((A).dev < (B).dev) : ((A).Lab[0] > (B).Lab[0]))
+ HEAPSORT(wval, pvals[j], n_pvals[j]);
+#undef HEAP_COMPARE
+
+ /* Check the maximum value looks OK */
+ if (n_pvals[j] < 5)
+ warning("Channel %d has only %d test patches",n_pvals[j]);
+ if (pvals[j][n_pvals[j]-1].dev < 0.99)
+ warning("Channel %d has max test patch value of %f",pvals[j][n_pvals[j]-1]);
+
+ if (verb > 1) {
+ printf("Chan %d has %d raw values:\n",j,n_pvals[j]);
+ for (i = 0; i < n_pvals[j]; i++) {
+ wp = &pvals[j][i];
+ printf(" %d: dev %f,XYZ %f %f %f, Lab %f %f %f\n",
+ i,wp->dev,wp->XYZ[0], wp->XYZ[1], wp->XYZ[2],
+ wp->Lab[0], wp->Lab[1], wp->Lab[2]);
+ }
+ }
+ }
+
+ if (sp2cie != NULL)
+ sp2cie->del(sp2cie);
+ }
+ icg->del(icg); /* Clean up */
+
+ /* Interpolate Lab using rspl */
+ for (j = 0; j < devchan; j++) {
+ datai low,high;
+ datao olow,ohigh;
+ int gres[MXDI];
+ double avgdev[MXDO];
+ cow *dpoints;
+
+ low[0] = 0.0;
+ high[0] = 1.0;
+ gres[0] = GRES;
+ olow[0] = 0.0;
+ ohigh[0] = 100.0;
+ olow[1] = olow[2] = -128.0;
+ ohigh[1] = ohigh[2] = 128.0;
+ avgdev[0] = 0.0025;
+ avgdev[1] = 0.005;
+ avgdev[2] = 0.005;
+
+ if ((raw[j] = new_rspl(RSPL_NOFLAGS,1, 3)) == NULL)
+ error("new_rspl() failed");
+
+ if ((dpoints = (cow *)malloc(sizeof(cow) * n_pvals[j])) == NULL)
+ error("malloc dpoints[%d] failed",n_pvals[j]);
+
+ for (i = 0; i < n_pvals[j]; i++) {
+ dpoints[i].p[0] = pvals[j][i].dev;
+ dpoints[i].v[0] = pvals[j][i].Lab[0];
+ dpoints[i].v[1] = pvals[j][i].Lab[1];
+ dpoints[i].v[2] = pvals[j][i].Lab[2];
+ if (i == 0)
+ dpoints[i].w = (double)n_white;
+ else
+ dpoints[i].w = 1.0;
+ }
+
+ raw[j]->fit_rspl_w(raw[j],
+ RSPLFLAGS,
+ dpoints, /* Test points */
+ n_pvals[j], /* Number of test points */
+ low, high, gres, /* Low, high, resolution of grid */
+ olow, ohigh, /* Default data scale */
+ smooth, /* Smoothing */
+ avgdev, /* Average deviation */
+ NULL); /* iwidth */
+
+
+ if (verb > 1) {
+ }
+
+ free(dpoints);
+ }
+
+ /* Plot the raw curves */
+ if (doplot > 1) {
+ double xx[PRES];
+ double yy[3][PRES];
+
+ for (j = 0; j < devchan; j++) {
+ printf("Chan %d raw L*a*b*:\n",j);
+ for (i = 0; i < PRES; i++) {
+ co tp; /* Test point */
+ xx[i] = i/(double)(PRES-1);
+ tp.p[0] = xx[i];
+ raw[j]->interp(raw[j], &tp);
+ yy[0][i] = tp.v[0];
+ yy[1][i] = tp.v[1];
+ yy[2][i] = tp.v[2];
+ }
+ do_plot(xx, yy[0], yy[1], yy[2], PRES);
+ }
+ }
+
+ /* Create a RSPL of absolute deltaE and relative deltaE '94 */
+ for (j = 0; j < devchan; j++) {
+ datai low,high;
+ int gres[MXDI];
+ double avgdev[MXDO];
+ co *dpoints_a;
+ co *dpoints_r;
+ double wh[3], prev[3], tot;
+
+ low[0] = 0.0;
+ high[0] = 1.0;
+ gres[0] = GRES;
+ avgdev[0] = 0.0;
+
+ if ((ade[j] = new_rspl(RSPL_NOFLAGS,1, 1)) == NULL)
+ error("new_rspl() failed");
+ if (imitate) {
+ if ((pcade[j] = new_rspl(RSPL_NOFLAGS,1, 1)) == NULL)
+ error("new_rspl() failed");
+ }
+ if ((rde[j] = new_rspl(RSPL_NOFLAGS,1, 1)) == NULL)
+ error("new_rspl() failed");
+
+ if ((dpoints_a = malloc(sizeof(co) * GRES)) == NULL)
+ error("malloc dpoints[%d] failed",GRES);
+ if ((dpoints_r = malloc(sizeof(co) * GRES)) == NULL)
+ error("malloc dpoints[%d] failed",GRES);
+
+//printf("~1 Chan %d:\n",j);
+ for (i = 0; i < GRES; i++) {
+ co tp; /* Test point */
+
+ tp.p[0] = i/(double)(GRES-1);
+ raw[j]->interp(raw[j], &tp);
+
+ dpoints_a[i].p[0] = tp.p[0];
+ dpoints_r[i].p[0] = tp.p[0];
+ if (i == 0) {
+//printf("~1 wht = %f %f %f\n",tp.v[0],tp.v[1],tp.v[2]);
+ tot = 0.0;
+ prev[0] = wh[0] = tp.v[0];
+ prev[1] = wh[1] = tp.v[1];
+ prev[2] = wh[2] = tp.v[2];
+ dpoints_a[i].v[0] = 0.0;
+ dpoints_r[i].v[0] = 0.0;
+ } else {
+//printf("~1 samp %d = %f %f %f\n",i,tp.v[0],tp.v[1],tp.v[2]);
+ /* Use Euclidean for large DE: (CIE94 stuffs up here) */
+ dpoints_a[i].v[0] = icmLabDE(tp.v, wh);
+ /* And CIE94 for small: */
+ tot += icmCIE94(tp.v, prev);
+ prev[0] = tp.v[0];
+ prev[1] = tp.v[1];
+ prev[2] = tp.v[2];
+ dpoints_r[i].v[0] = tot;
+ }
+//printf("~1 %d: dev %f, ade %f, rde %f\n",i,tp.p[0],dpoints_a[i].v[0],dpoints_r[i].v[0]);
+ }
+
+ ade[j]->set_rspl(ade[j],
+ 0,
+ (void *)dpoints_a, /* Test points */
+ rsplset1, /* Setting function */
+ low, high, gres, /* Low, high, resolution of grid */
+ NULL, NULL /* Default data scale */
+ );
+ if (imitate) {
+ pcade[j]->set_rspl(pcade[j],
+ 0,
+ (void *)dpoints_a, /* Test points */
+ rsplset1, /* Setting function */
+ low, high, gres, /* Low, high, resolution of grid */
+ NULL, NULL /* Default data scale */
+ );
+ }
+ rde[j]->set_rspl(rde[j],
+ 0,
+ (void *)dpoints_r, /* Test points */
+ rsplset1, /* Setting function */
+ low, high, gres, /* Low, high, resolution of grid */
+ NULL, NULL /* Default data scale */
+ );
+ free(dpoints_a);
+ free(dpoints_r);
+ }
+
+ if (initial) {
+ /* Establish the ademax values */
+ pct->update_devmax(pct, -1, -1.0); /* Make sure there is a value for each */
+ pct->update_ademax(pct, -1, -1.0);
+ pct->update_ademin(pct, -1, -1.0);
+ for (j = 0; j < devchan; j++) {
+ co tp; /* Test point */
+
+ if (pct->devmax[j] < 0.0) { /* Auto */
+ double maxd, maxde, maxix;
+ /* Locate the point of maximum aDE */
+ for (maxde = -1.0, i = 0; i < GRES; i++) {
+
+ tp.p[0] = i/(GRES-1.0);
+ ade[j]->interp(ade[j], &tp);
+ if (tp.v[0] > maxde) {
+ maxd = tp.p[0];
+ maxde = tp.v[0];
+ maxix = i;
+ }
+ }
+ pct->devmax[j] = maxd;
+ pct->ademax[j] = maxde; /* Temporary */
+//printf("Chan %d, dev %f, max de = %f\n", j, maxd, maxde);
+
+ if (maxd < 0.2) {
+ warning("Chan %d, max DE point %f is below < 0.2 - ignored\n", j, maxd);
+ maxix = GRES-1;
+ }
+ /* Then locate the point below that where the slope */
+ /* becomes reasonable. */
+ for (i = maxix; i >= 40; i--) {
+ double aslope, minslope = 1e6;
+ double naslope, nminslope;
+ int k;
+
+ /* Compute the minimum over a span of 20/GRES */
+ for (k = 0; k < 40; k++) {
+ double dp, dv, slope;
+
+ tp.p[0] = (i-k)/(GRES-1.0);
+ dp = tp.p[0];
+ ade[j]->interp(ade[j], &tp);
+ dv = tp.v[0];
+
+ tp.p[0] = (i-k-1)/(GRES-1.0);
+ ade[j]->interp(ade[j], &tp);
+ slope = (dv - tp.v[0])/(dp - (i-k-1)/(GRES-1.0));
+ if (k == 0)
+ aslope = slope;
+//printf(" Chan %d, dev %f, dv = %f, slope = %f\n", j, (i-k)/(GRES-1.0),dv - tp.v[0],slope);
+ if (slope < minslope)
+ minslope = slope;
+ }
+//printf("Chan %d, dev %f, aslope = %f, min slope = %f\n", j, i/(GRES-1.0),aslope,minslope);
+ /* Normalize the slopes */
+ naslope = aslope * SLOPE_NORM/pct->ademax[j];
+ nminslope = minslope * SLOPE_NORM/pct->ademax[j];
+//printf("Chan %d, dev %f, norm aslope = %f, min slope = %f\n", j, i/(GRES-1.0),naslope,nminslope);
+
+ if (naslope > MIN_SLOPE_A && nminslope >= MIN_SLOPE_O)
+ break;
+ }
+ pct->devmax[j] = i/(GRES-1.0);
+
+ /* Scale auto max device value */
+ if (maxscale[j] >= 0.0)
+ pct->devmax[j] *= maxscale[j];
+
+ /* Manually set initial dev max */
+ } else {
+ if (maxscale[j] >= 0.0)
+ warning("Chan %d, scale %.1f%% of auto max ignored since max override used\n", j, maxscale[j] * 100.0);
+ }
+
+ /* Lookup devmax to set ademax */
+ tp.p[0] = pct->devmax[j];
+ ade[j]->interp(ade[j], &tp);
+ pct->ademax[j] = tp.v[0];
+
+
+ /* Establish a default ademin value */
+ if (pct->ademin[j] < 0.0)
+ pct->ademin[j] = 0.0;
+ }
+
+ } else if (recal) {
+
+ /* Since the plot markers use devmax, look it up */
+ for (j = 0; j < devchan; j++) {
+ if ((pct->devmax[j] = rspl_ilookup(ade[j], 0.5, pct->ademax[j])) < 0.0)
+ error("Unexpected failure to invert curve %d for ADE %f",j,pct->ademax[j]);
+ }
+ }
+
+ /* Find the maximum aDE value for each curve */
+ for (j = 0; j < devchan; j++) {
+ co tp; /* Test point */
+ mxade[j] = -1e6;
+ for (i = 0; i < PRES; i++) {
+ tp.p[0] = i/(double)(PRES-1);
+ ade[j]->interp(ade[j], &tp);
+ if (tp.v[0] > mxade[j])
+ mxade[j] = tp.v[0];
+ }
+ }
+
+ if (initial || recal) {
+ /* Compute an ideal power-like value for test target */
+ for (j = 0; j < devchan; j++) {
+ double hdv; /* Half device value */
+ double hdvrde; /* Half device rDE value */
+ double thdvrde; /* Target half device rDE value */
+ double thdv; /* Target half device value */
+ double fdvrde; /* Full device rDE value */
+ co tp; /* Test point */
+
+ /* full rDE */
+ tp.p[0] = pct->devmax[j];
+ rde[j]->interp(rde[j], &tp);
+ fdvrde = tp.v[0];
+
+ /* Half device value of maximum */
+ hdv = pct->devmax[j] * 0.5;
+
+ /* rDE value half the device value */
+ tp.p[0] = hdv;
+ rde[j]->interp(rde[j], &tp);
+ hdvrde = tp.v[0];
+
+ /* rDE value we'd like at half the device value */
+ thdvrde = 0.5 * fdvrde;
+
+ /* Device value to get the rDE value we'd like at half */
+ if ((thdv = rspl_ilookup(rde[j], 0.5, thdvrde)) < 0.0)
+ error("Unexpected failure to invert curve %d for ADE %f",j,thdvrde);
+
+//printf("hdv %f, hdvrde %f, thdvrde %f, fdvrde %f, thdv %f\n",hdv,hdvrde,thdvrde, fdvrde,thdv);
+
+ /* Power like value needed to get rDE value we'd like at hald device */
+ idpow[j] = icx_powlike_needed(hdv, thdv);
+
+ }
+ }
+
+ if (verb > 1) {
+ printf("Abs DE values:\n");
+ for (i = 0; i < PRES; i++) {
+ co tp; /* Test point */
+ tp.p[0]= i/(double)(PRES-1);
+
+ printf(" dev %f, aDE",tp.p[0]);
+ for (j = 0; j < 6 && j < devchan; j++) {
+ ade[j]->interp(ade[j], &tp);
+ printf(" %f",tp.v[0]);
+ }
+ printf("\n");
+ }
+ printf("Rel DE values:\n");
+ for (i = 0; i < PRES; i++) {
+ co tp; /* Test point */
+ tp.p[0]= i/(double)(PRES-1);
+
+ printf(" dev %f, rdev",tp.p[0]);
+ for (j = 0; j < 6 && j < devchan; j++) {
+ rde[j]->interp(rde[j], &tp);
+ printf(" %f",tp.v[0]);
+ }
+ printf("\n");
+ }
+ }
+
+ if (initial || recal) {
+ if (verb && pct->is_set(pct)) {
+ for (j = 0; j < devchan; j++) {
+ printf("Chan %d Dev max %f, aDE Max %f, aDE Min %f\n",j,pct->devmax[j],pct->ademax[j],pct->ademin[j]);
+ }
+ }
+
+ if (verb) {
+ double avgpow = 0.0;
+ for (j = 0; j < devchan; j++) {
+ printf("Chan %d ideal targen power = %f\n",j,idpow[j]);
+ avgpow += idpow[j];
+ }
+ avgpow /= (double)devchan;
+ printf("Average ideal targen power = %f\n",avgpow);
+ }
+ }
+
+ /* Plot both the delta E curves, and markers */
+ if (doplot) {
+ co tp; /* Test point */
+ double xx[PRES];
+ double yy[10][PRES];
+ double cx[10], cy[10];
+ int nmark;
+
+ printf("Absolute DE plot:\n");
+
+ for (i = 0; i < PRES; i++) {
+ xx[i] = i/(double)(PRES-1);
+
+ for (j = 0; j < 10 && j < devchan; j++) {
+ tp.p[0] = xx[i];
+ ade[j]->interp(ade[j], &tp);
+ yy[j][i] = tp.v[0];
+ }
+ }
+ nmark = 0;
+ if (pct->is_set(pct)) {
+ /* Add markers for deMax */
+ for (j = 0; j < 10 && j < devchan; j++) {
+ cx[j] = pct->devmax[j];
+ cy[j] = pct->ademax[j];
+ nmark++;
+ }
+ }
+ do_plot10p(xx, devchan > 3 ? yy[3] : NULL,
+ devchan > 1 ? yy[1] : NULL,
+ devchan > 4 ? yy[4] : NULL,
+ devchan > 0 ? yy[0] : NULL,
+ devchan > 2 ? yy[2] : NULL,
+ devchan > 5 ? yy[5] : NULL,
+ devchan > 6 ? yy[6] : NULL,
+ devchan > 7 ? yy[7] : NULL,
+ devchan > 8 ? yy[8] : NULL,
+ devchan > 9 ? yy[9] : NULL,
+ PRES,
+ cx, cy, verify ? 0 : nmark);
+
+ printf("Relative DE plot:\n");
+
+ for (i = 0; i < PRES; i++) {
+ xx[i] = i/(double)(PRES-1.0);
+
+ for (j = 0; j < 10 && j < devchan; j++) {
+ tp.p[0] = xx[i];
+ rde[j]->interp(rde[j], &tp);
+ yy[j][i] = tp.v[0];
+ }
+ }
+ nmark = 0;
+ if (pct->is_set(pct)) {
+ /* Add markers for deMax */
+ for (j = 0; j < 10 && j < devchan; j++) {
+ cx[j] = pct->devmax[j];
+ tp.p[0] = cx[j];
+ rde[j]->interp(rde[j], &tp);
+ cy[j] = tp.v[0];
+ nmark++;
+ }
+ }
+ do_plot10p(xx, devchan > 3 ? yy[3] : NULL,
+ devchan > 1 ? yy[1] : NULL,
+ devchan > 4 ? yy[4] : NULL,
+ devchan > 0 ? yy[0] : NULL,
+ devchan > 2 ? yy[2] : NULL,
+ devchan > 5 ? yy[5] : NULL,
+ devchan > 6 ? yy[6] : NULL,
+ devchan > 7 ? yy[7] : NULL,
+ devchan > 8 ? yy[8] : NULL,
+ devchan > 9 ? yy[9] : NULL,
+ PRES,
+ cx, cy, verify ? 0 : nmark);
+
+ if (idpow[0] > 0.0) {
+ printf("Relative DE plot with ideal targen power applied:\n");
+
+ for (i = 0; i < PRES; i++) {
+ xx[i] = i/(double)(PRES-1.0);
+
+ for (j = 0; j < 10 && j < devchan; j++) {
+ tp.p[0] = icx_powlike(xx[i],idpow[j]);
+ rde[j]->interp(rde[j], &tp);
+ yy[j][i] = tp.v[0];
+ }
+ }
+ nmark = 0;
+ if (pct->is_set(pct)) {
+ /* Add markers for deMax */
+ for (j = 0; j < 10 && j < devchan; j++) {
+ cx[j] = pct->devmax[j];
+ tp.p[0] = icx_powlike(cx[j],idpow[j]);
+ rde[j]->interp(rde[j], &tp);
+ cy[j] = tp.v[0];
+ nmark++;
+ }
+ }
+ do_plot10p(xx, devchan > 3 ? yy[3] : NULL,
+ devchan > 1 ? yy[1] : NULL,
+ devchan > 4 ? yy[4] : NULL,
+ devchan > 0 ? yy[0] : NULL,
+ devchan > 2 ? yy[2] : NULL,
+ devchan > 5 ? yy[5] : NULL,
+ devchan > 6 ? yy[6] : NULL,
+ devchan > 7 ? yy[7] : NULL,
+ devchan > 8 ? yy[8] : NULL,
+ devchan > 9 ? yy[9] : NULL,
+ PRES,
+ cx, cy, verify ? 0 : nmark);
+ }
+ }
+
+ /* Compare the previous expected aDE against the current one */
+ if (verify) {
+ co tp; /* Test point */
+ double avg[MAX_CHAN];
+ double max[MAX_CHAN];
+ double rms[MAX_CHAN];
+ int verified = 1;
+
+ /* Verify each channel */
+ for (j = 0; j < devchan; j++) {
+ co tp;
+
+ avg[j] = 0.0;
+ max[j] = -1.0;
+ rms[j] = 0.0;
+
+ /* Sample it at GRES */
+ for (i = 0; i < GRES; i++) {
+ double iv, targ, val, tt;
+
+ iv = i/(GRES-1.0);
+
+ /* Lookup the ade that we expect */
+ tp.p[0] = iv;
+ pcade[j]->interp(pcade[j], &tp);
+ targ = tp.v[0];
+
+ /* Lookup the ade that we have */
+ tp.p[0] = iv;
+ ade[j]->interp(ade[j], &tp);
+ val = tp.v[0];
+
+ /* Compute the stats */
+ tt = fabs(targ - val);
+ avg[j] += tt;
+ rms[j] += tt * tt;
+ if (tt > max[j])
+ max[j] = tt;
+//printf("~1 chan %d, ix %d, iv %f, targ %f, actual %f, err %f\n",j,i,iv,targ,val,tt);
+ }
+ avg[j] /= (double)GRES;
+ rms[j] /= (double)GRES;
+ rms[j] = sqrt(rms[j]);
+ if (max[j] > ver_maxde)
+ verified = 0;
+ }
+ if (verb) {
+ for (j = 0; j < devchan; j++) {
+ printf("Verify results:\n");
+ printf("Channel %d has DE avg %.1f, rms %.1f, max %.1f\n",j,avg[j],rms[j],max[j]);
+ }
+ if (verified)
+ printf("Verified OK\n");
+ else
+ printf("Verification FAILED\n");
+ }
+
+ /* Plot the verification curves */
+ if (doplot) {
+ double xx[PRES];
+ double yy[10][PRES];
+
+ printf("Verification match plot:\n");
+
+ for (j = 0; j < 6 && j < devchan; j++) {
+ co tp; /* Test point */
+ double max;
+
+ /* Establish the scale */
+ tp.p[0] = 1.0;
+ pcade[j]->interp(pcade[j], &tp);
+ max = tp.v[0];
+
+ for (i = 0; i < PRES; i++) {
+ xx[i] = i/(double)(PRES-1);
+
+ /* Convert ade target to device */
+ tp.v[0] = max * xx[i];
+ if ((tp.p[0] = rspl_ilookup(pcade[j], 1.0, tp.v[0])) < 0.0)
+ error("Unexpected failure to invert curve %d for pcADE %f",j,tp.v[0]);
+ /* Convert device to actual ade */
+ ade[j]->interp(ade[j], &tp);
+ if (fabs(max) > 0.1)
+ yy[j][i] = tp.v[0]/max;
+ else
+ yy[j][i] = 0.0;
+ }
+ }
+ do_plot10(xx, devchan > 3 ? yy[3] : NULL,
+ devchan > 1 ? yy[1] : NULL,
+ devchan > 4 ? yy[4] : NULL,
+ devchan > 0 ? yy[0] : NULL,
+ devchan > 2 ? yy[2] : NULL,
+ devchan > 5 ? yy[5] : NULL,
+ devchan > 6 ? yy[6] : NULL,
+ devchan > 7 ? yy[7] : NULL,
+ devchan > 8 ? yy[8] : NULL,
+ devchan > 9 ? yy[9] : NULL,
+ PRES, 0);
+ if (!verified)
+ exit(1);
+ }
+ } else if (initial) {
+
+ /* Convert any transfer curve target points into a smooth curve */
+ if (pct->no_tpoints > 0) {
+ int gres[MXDI] = { GRES };
+ co *pnts;
+
+ if ((pnts = (co *)calloc(pct->no_tpoints + 2, sizeof(co))) == NULL)
+ error ("Malloc of rspl points failed");
+
+ for (j = 0; j < devchan; j++) {
+ int npts;
+ int gotmin, gotmax;
+
+ /* Count the number of valid points */
+ for (npts = i = 0; i < pct->no_tpoints; i++) {
+ if (pct->tpoints[i].val[j] >= 0.0)
+ npts++;
+ }
+ if (npts == 0)
+ continue; /* No target curve for this channel */
+
+ if ((tcurves[j] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL)
+ error("new_rspl(1,1) failed");
+
+ gotmin = gotmax = 0;
+ for (npts = i = 0; i < pct->no_tpoints; i++) {
+ if (pct->tpoints[i].val[j] < 0.0)
+ continue;
+
+ pnts[npts].p[0] = pct->tpoints[i].loc;
+ pnts[npts].v[0] = pct->tpoints[i].val[j];
+ if (pnts[npts].p[0] < 0.0)
+ pnts[npts].p[0] = 0.0;
+ else if (pnts[npts].p[0] > 1.0)
+ pnts[npts].p[0] = 1.0;
+ if (pnts[npts].v[0] < 0.0)
+ pnts[npts].v[0] = 0.0;
+ else if (pnts[npts].v[0] > 1.0)
+ pnts[npts].v[0] = 1.0;
+
+ if (pnts[npts].p[0] <= 0.05)
+ gotmin = 1;
+ if (pnts[npts].p[0] >= 0.95)
+ gotmax = 1;
+ npts++;
+ }
+ /* Add default anchors if there are none supplied */
+ if (gotmin == 0) {
+ pnts[npts].p[0] = 0.0;
+ pnts[npts++].v[0] = 0.0;
+ }
+ if (gotmax == 0) {
+ pnts[npts].p[0] = 1.0;
+ pnts[npts++].v[0] = 1.0;
+ }
+
+ /* Fit the curve to the given points */
+ tcurves[j]->fit_rspl(tcurves[j], RSPLFLAGS, pnts, npts,
+ NULL, NULL, gres, NULL, NULL, TCURVESMOOTH, NULL, NULL);
+ }
+ free(pnts);
+
+ /* Plot the target curves */
+ if (doplot) {
+ double xx[PRES];
+ double yy[10][PRES];
+
+ printf("Target curves plot:\n");
+
+ for (i = 0; i < (PRES-1); i++) {
+ co tp; /* Test point */
+ double pp = i/(PRES-1.0);
+
+ xx[i] = pp;
+ for (j = 0; j < 10 && j < devchan; j++) {
+ if (tcurves[j] != NULL) {
+ tp.p[0] = pp;
+ tcurves[j]->interp(tcurves[j], &tp);
+ yy[j][i] = tp.v[0];
+ } else
+ yy[j][i] = pp;
+ }
+ }
+ do_plot10(xx, devchan > 3 ? yy[3] : NULL,
+ devchan > 1 ? yy[1] : NULL,
+ devchan > 4 ? yy[4] : NULL,
+ devchan > 0 ? yy[0] : NULL,
+ devchan > 2 ? yy[2] : NULL,
+ devchan > 5 ? yy[5] : NULL,
+ devchan > 6 ? yy[6] : NULL,
+ devchan > 7 ? yy[7] : NULL,
+ devchan > 8 ? yy[8] : NULL,
+ devchan > 9 ? yy[9] : NULL,
+ PRES, 0);
+ }
+ }
+
+ /* Do inverse lookup to create relative linearization curves */
+ n_cvals = CAL_RES;
+ for (j = 0; j < devchan; j++) {
+ co tp;
+ double rdemin, rdemax; /* Relative DE min and max targets */
+
+ /* Convert absolute de aims to relative */
+ if ((rdemax = rspl_ilookup(ade[j], 0.0, pct->ademax[j])) < 0.0)
+ error("Unexpected failure to invert curve %d for DE %f",j,pct->ademax[j]);
+
+ tp.p[0] = rdemax;
+ rde[j]->interp(rde[j], &tp);
+ rdemax = tp.v[0];
+
+ if ((rdemin = rspl_ilookup(ade[j], 1.0, pct->ademin[j])) < 0.0)
+ error("Unexpected failure to invert curve %d for DE %f",j,pct->ademax[j]);
+
+ tp.p[0] = rdemin;
+ rde[j]->interp(rde[j], &tp);
+ rdemin = tp.v[0];
+
+ if (verb > 0)
+ printf("Chan %d: rDE Max = %f, rDE Min = %f\n",j,rdemax,rdemin);
+
+ if ((cvals[j] = (wval *)malloc(sizeof(wval) * n_cvals)) == NULL)
+ error("Malloc of %d cvals failed",n_cvals);
+
+ /* Convert relative delta E aim to device value */
+ for (i = 0; i < n_cvals; i++) {
+ double x = i/(n_cvals-1.0);
+ double inv;
+
+ cvals[j][i].inv = x;
+
+ /* Apply any aim tweak curve */
+ if (tcurves[j] != NULL) {
+ tp.p[0] = x;
+ tcurves[j]->interp(tcurves[j], &tp);
+ x = tp.v[0];
+ }
+
+ inv = x * (rdemax - rdemin) + rdemin;
+
+ if ((cvals[j][i].dev = rspl_ilookup(rde[j], 0.5, inv)) < 0.0)
+ error("Unexpected failure to invert curve %d for DE %f",j,inv);
+//printf("~1 chan %d, step %d, inv %f, detarg %f, got dev %f\n",j,i,x,pp[0].v[0],pp[k].p[0]);
+ }
+ }
+
+ } else if (recal || imitate) {
+
+ n_cvals = CAL_RES;
+ for (j = 0; j < devchan; j++) {
+ co tp;
+
+ if ((cvals[j] = (wval *)malloc(sizeof(wval) * n_cvals)) == NULL)
+ error("Malloc of %d cvals failed",n_cvals);
+
+ /* Lookup the expected ade for each input device value, and */
+ /* then translate it into the required output device value */
+ for (i = 0; i < n_cvals; i++) {
+ double x = i/(n_cvals-1.0);
+
+ cvals[j][i].inv = tp.p[0] = x;
+ pcade[j]->interp(pcade[j], &tp);
+
+ if ((cvals[j][i].dev = rspl_ilookup(ade[j], 0.5, tp.v[0])) < 0.0)
+ error("Unexpected failure to invert curve %d for DE %f",j,tp.v[0]);
+//printf("~1 chan %d, ix %d, inv %f, pcade %f, iade %f\n",j,i,x,tp.v[0],cvals[j][i].dev);
+ }
+ }
+ }
+
+ if (initial || recal || imitate) {
+
+ if (verb > 1) {
+ printf("Calibration curve values:\n");
+ for (i = 0; i < n_cvals; i++) {
+ printf(" inv %f, dev",cvals[0][i].inv);
+ for (j = 0; j < devchan; j++) {
+ printf(" %f",cvals[j][i].dev);
+ }
+ printf("\n");
+ }
+ }
+
+ /* Plot the calibration curves */
+ if (doplot) {
+ double xx[PRES];
+ double yy[10][PRES];
+
+ printf("Calibration curve plot:\n");
+
+ for (i = 0; i < n_cvals; i++) {
+ xx[i] = cvals[0][i].inv;
+ for (j = 0; j < 10 && j < devchan; j++) {
+ yy[j][i] = cvals[j][i].dev;
+ }
+ }
+ do_plot10(xx, devchan > 3 ? yy[3] : NULL, /* Black */
+ devchan > 1 ? yy[1] : NULL, /* Red */
+ devchan > 4 ? yy[4] : NULL, /* Green */
+ devchan > 0 ? yy[0] : NULL, /* Blue */
+ devchan > 2 ? yy[2] : NULL, /* Yellow */
+ devchan > 5 ? yy[5] : NULL, /* Purple */
+ devchan > 6 ? yy[6] : NULL, /* Brown */
+ devchan > 7 ? yy[7] : NULL, /* Orange */
+ devchan > 8 ? yy[8] : NULL, /* Grey */
+ devchan > 9 ? yy[9] : NULL, /* White */
+ PRES, 0);
+ }
+
+ /* Write out an Argyll .CAL file */
+ if (dowrite) {
+ cgats *ocg; /* output cgats structure */
+ time_t clk = time(0);
+ struct tm *tsp = localtime(&clk);
+ char *atm = asctime(tsp); /* Ascii time */
+ char *ident = icx_inkmask2char(devmask, 1);
+ char *bident = icx_inkmask2char(devmask, 0);
+ cgats_set_elem *setel; /* Array of set value elements */
+ int nsetel = 0;
+ int ncps; /* Number of curve parameters */
+ double *cps[3]; /* Arrays of curve parameters */
+ char *bp = NULL, buf[100]; /* Buffer to sprintf into */
+ co tp;
+
+ ocg = new_cgats(); /* Create a CGATS structure */
+ ocg->add_other(ocg, "CAL"); /* our special type is Calibration file */
+
+ ocg->add_table(ocg, tt_other, 0); /* Add a table for RAMDAC values */
+ ocg->add_kword(ocg, 0, "DESCRIPTOR", "Argyll Device Calibration Curves",NULL);
+ ocg->add_kword(ocg, 0, "ORIGINATOR", "Argyll printcal", NULL);
+ atm[strlen(atm)-1] = '\000'; /* Remove \n from end */
+ ocg->add_kword(ocg, 0, "CREATED",atm, NULL);
+
+ ocg->add_kword(ocg, 0, "DEVICE_CLASS","OUTPUT", NULL);
+ ocg->add_kword(ocg, 0, "COLOR_REP", ident, NULL);
+
+ if (xpi.deviceMfgDesc != NULL)
+ ocg->add_kword(ocg, 0, "MANUFACTURER", xpi.deviceMfgDesc, NULL);
+ if (xpi.modelDesc != NULL)
+ ocg->add_kword(ocg, 0, "MODEL", xpi.modelDesc, NULL);
+ if (xpi.profDesc != NULL)
+ ocg->add_kword(ocg, 0, "DESCRIPTION", xpi.profDesc, NULL);
+ if (xpi.copyright != NULL)
+ ocg->add_kword(ocg, 0, "COPYRIGHT", xpi.copyright, NULL);
+
+ /* Setup the table which holds the translation from calibrated */
+ /* device value "I" to the raw device channel value */
+ sprintf(buf, "%s_I",bident);
+ ocg->add_field(ocg, 0, buf, r_t);
+ nsetel++;
+ for (j = 0; j < devchan; j++) {
+ inkmask imask = icx_index2ink(devmask, j);
+ sprintf(buf, "%s_%s",bident,icx_ink2char(imask));
+ ocg->add_field(ocg, 0, buf, r_t);
+ nsetel++;
+ }
+ if ((setel = (cgats_set_elem *)malloc(sizeof(cgats_set_elem) * nsetel)) == NULL)
+ error("Malloc failed!");
+
+ /* Write the per channel device to device loolup curve values */
+ if (devmask & ICX_ADDITIVE) {
+ for (i = n_cvals-1; i >= 0; i--) {
+
+ setel[0].d = 1.0 - cvals[0][i].inv;
+ for (j = 0; j < devchan; j++)
+ setel[1+j].d = 1.0 - cvals[j][i].dev;
+ ocg->add_setarr(ocg, 0, setel);
+ }
+ } else {
+ for (i = 0; i < n_cvals; i++) {
+
+ setel[0].d = cvals[0][i].inv;
+ for (j = 0; j < devchan; j++)
+ setel[1+j].d = cvals[j][i].dev;
+ ocg->add_setarr(ocg, 0, setel);
+ }
+ }
+
+ free(setel);
+
+ /* Write the calibration target information to a second table */
+ if (pct->write(pct, ocg, 1) != 0)
+ error("Writing cal target info to cal file '%s'",outname);
+
+ /* Add a third table which is the expected absolute DE response */
+ /* of the calibrated device. */
+ ocg->add_table(ocg, tt_other, 0); /* Add a table for RAMDAC values */
+ ocg->add_kword(ocg, 2, "DESCRIPTOR", "Argyll Output Calibration Expected DE Response",NULL);
+ ocg->add_kword(ocg, 2, "ORIGINATOR", "Argyll printcal", NULL);
+ atm[strlen(atm)-1] = '\000'; /* Remove \n from end */
+ ocg->add_kword(ocg, 2, "CREATED",atm, NULL);
+
+ ocg->add_kword(ocg, 2, "DEVICE_CLASS","OUTPUT", NULL);
+ ocg->add_kword(ocg, 2, "COLOR_REP", ident, NULL);
+
+ /* Setup the table which holds the translation from calibrated */
+ /* device value "I" to the expected absolute deltaE value for */
+ /* each colorant. */
+ sprintf(buf, "%s_I_DE",bident);
+ ocg->add_field(ocg, 2, buf, r_t);
+ nsetel++;
+ for (j = 0; j < devchan; j++) {
+ inkmask imask = icx_index2ink(devmask, j);
+ sprintf(buf, "%s_%s_DE",bident,icx_ink2char(imask));
+ ocg->add_field(ocg, 2, buf, r_t);
+ nsetel++;
+ }
+ if ((setel = (cgats_set_elem *)malloc(sizeof(cgats_set_elem) * nsetel)) == NULL)
+ error("Malloc failed!");
+
+ for (i = 0; i < n_cvals; i++) {
+ int ix = i;
+
+ /* Calibrated device value */
+ if (devmask & ICX_ADDITIVE) {
+ ix = n_cvals -1 - i;
+ setel[0].d = 1.0 - cvals[0][ix].inv;
+ } else
+ setel[0].d = cvals[0][ix].inv;
+
+ for (j = 0; j < devchan; j++) {
+ tp.p[0] = cvals[j][ix].dev; /* Raw device value */
+ ade[j]->interp(ade[j], &tp); /* Corresponding ade value */
+ setel[1+j].d = tp.v[0];
+ }
+ ocg->add_setarr(ocg, 2, setel);
+ }
+
+ free(setel);
+
+ if (ocg->write_name(ocg, outname))
+ error("Write error to file '%s': %s",outname,ocg->err);
+
+ if (verb)
+ printf("Written calibration file '%s'\n",outname);
+
+ ocg->del(ocg); /* Clean up */
+ free(ident);
+ free(bident);
+ }
+
+ /*
+ The structure of *.AMP is very simple.
+
+ It has 5 tables that have
+ 256 entries that are 8-bit (byte), even for 16 bit mode.
+
+ The first table 0000..00FFh is the CMYK channel.
+ The second table 0100..01FFh is the C channel.
+ The third table 0200..02FFh is the M channel.
+ The fourth table 0300..03FFh is the Y channel.
+ The fifth table 0300..03FFh is the K channel.
+
+ Table position nn00h is the black end and table position nnFFh
+ is the white end.
+
+ */
+
+ /* Write an Adobe map format (.AMP) file */
+ /* (It's not clear if more than 4 channels is allowed) */
+ if (dowrite && doamp) {
+ FILE *fp;
+ cgatsFile *p;
+ char nmode[50] = { '\000' };
+
+ strcpy(nmode, "w");
+#if defined(O_BINARY) || defined(_O_BINARY)
+ strcat(nmode, "b");
+#endif
+ if ((fp = fopen(ampname, nmode)) == NULL)
+ error("Couldn't open '%s' for writing",ampname);
+
+ /* CMYK table is unity */
+ for (i = 0; i < 256; i++) {
+ if (putc(i,fp) == EOF)
+ error("Error writing to fle '%s'",ampname);
+ }
+ for (j = 0; j < devchan; j++) {
+ for (i = 0; i < 256; i++) {
+ int x;
+ if (devmask & ICX_ADDITIVE)
+ x = (int)(cvals[j][i].dev * 255.0 + 0.5); /* ??? */
+ else
+ x = 255 - (int)(cvals[j][255 - i].dev * 255.0 + 0.5);
+ if (putc(x,fp) == EOF)
+ error("Error writing to fle '%s'",ampname);
+//printf("~1 chan %d, inv %d, dev %d\n",j,i,x);
+ }
+ }
+ /* Extra 1:1 table */
+ for (i = 0; i < 256; i++) {
+ if (putc(i,fp) == EOF)
+ error("Error writing to fle '%s'",ampname);
+ }
+
+ if (fclose(fp) != 0)
+ error("Closing '%s' failed",ampname);
+
+ if (verb)
+ printf("Written calibration curves to '%s'\n",ampname);
+ }
+ }
+
+ /* Free up various possible allocations */
+ for (j = 0; j < devchan; j++) {
+ if (pvals[j] != NULL)
+ free(pvals[j]);
+ if (raw[j] != NULL)
+ raw[j]->del(raw[j]);
+ if (ade[j] != NULL)
+ ade[j]->del(ade[j]);
+ if (rde[j] != NULL)
+ rde[j]->del(rde[j]);
+ if (pcade[j] != NULL)
+ pcade[j]->del(pcade[j]);
+ if (cvals[j] != NULL)
+ free(cvals[j]);
+ if (tcurves[j] != NULL)
+ tcurves[j]->del(tcurves[j]);
+ }
+
+ return 0;
+}
+
+
+
+
+
+
+
+
+
+