Initial import of argyll version 1.5.1-8debian/1.5.1-8

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;
+}
+
+
+
+
+
+
+
+
+
+