state from 2014-09-07 13:30 MESZ

1 files changed, 658 insertions, 89 deletions

diff --git a/target/targen.c b/target/targen.c
index 09dcc35..c868dc4 100644
--- a/target/targen.c
+++ b/target/targen.c
@@ -27,6 +27,7 @@
 	then suppliment the measured patches. Would have to add another
 	set of measurement columns to .ti1 & .ti2 to carry the
 	already measured values through, or do clumbsy post merge ? 
+	(Latter is easiest).
 
 	Would be nice to be able to generate secondary
 	color ramps (ie. CMY for RGB space, RGB for CMYK space.)
@@ -103,11 +104,19 @@
 #define VRML_DIAG		/* Enable option to dump a VRML of the resulting full spread points */
 #undef ADDRECCLIPPOINTS	/* Add ink limited clipping points to regular grid */
 #define EMPH_NEUTRAL	/* Emphasise neutral axis, like CIE94 does */
-#define NEMPH_DEFAULT 0.5	/* Default emphasis == 2 x CIE94 */
+#define NEMPH_DEFAULT 0.5	/* Default neutral axis emphasis == 2 x CIE94 */
+#define XPOW_DEFAULT 1.0	/* Default extra device power value = none */
+#define DEMPH_DEFAULT 1.0	/* Default dark region emphasis == none */
 #define DEFANGLE 0.3333	/* For simdlat and simplat */
 #define SIMDLAT_TYPE SIMDLAT_BCC	/* Simdlat geometry type */
 #define MATCH_TOLL 1e-3	/* Tollerance of device value to consider a patch a duplicate */
 
+/* Display rise and fall time delay model. This is CRT like */
+#define DISPLAY_RISE_TIME 0.03		/* Assumed rise time to 90% of target level */ 
+#define DISPLAY_FALL_TIME 0.12		/* Assumed fall time to 90% of target level */
+#define DISPLAY_SETTLE_AIM 0.01		/* Aim for 1% of true level */
+#define DISPLAY_ABS_AIM 0.0001		/* Aim for .01% of true absolute level */
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
@@ -165,6 +174,8 @@ struct _pcpt {
 
 	/* Tuning parameters */
 	double nemph;		/* neutral emphasis, 0.0 - 1.0. Default 0.35 for == CIE94 */
+	double idemph;		/* inv. dark emphasis, 1.0 - 4.0. Default 1.0 == none */
+	double ixpow;		/* inv. extra power Default 1.0 == none */
 
 	/* ICC profile based */
 	icmFile *fp;
@@ -200,10 +211,10 @@ pcpt_to_XYZ(pcpt *s, double *out, double *in) {
 
 	if (s->xmask == s->nmask) {
 		for (e = 0; e < s->di; e++)
-			inv[e] = in[e];
+			inv[e] = icx_powlike(in[e], s->ixpow);
 	} else {
 		for (e = 0; e < s->di; e++)
-			inv[e] = 1.0 - in[e];
+			inv[e] = 1.0 - icx_powlike(in[e], s->ixpow);
 	}
 	if (s->luo2 != NULL)
 		s->luo2->lookup(s->luo2, out, inv);
@@ -230,10 +241,10 @@ pcpt_to_rLab(pcpt *s, double *out, double *in) {
 
 	if (s->xmask == s->nmask) {
 		for (e = 0; e < s->di; e++)
-			inv[e] = in[e];
+			inv[e] = icx_powlike(in[e], s->ixpow);
 	} else {
 		for (e = 0; e < s->di; e++)
-			inv[e] = 1.0 - in[e];
+			inv[e] = 1.0 - icx_powlike(in[e], s->ixpow);
 	}
 	if (s->luo != NULL)
 		s->luo->lookup(s->luo, out, inv);
@@ -251,6 +262,7 @@ pcpt_to_rLab(pcpt *s, double *out, double *in) {
 
 /* Perceptual conversion function */
 /* Internal device values 0.0 - 1.0 are converted into perceptually uniform 0.0 - 100.0 */
+/* This is used by optimal spread functions ? */
 static void
 pcpt_to_nLab(pcpt *s, double *out, double *in) {
 	int e;
@@ -258,10 +270,10 @@ pcpt_to_nLab(pcpt *s, double *out, double *in) {
 
 	if (s->xmask == s->nmask) {
 		for (e = 0; e < s->di; e++)
-			inv[e] = in[e];
+			inv[e] = icx_powlike(in[e], s->ixpow);
 	} else {
 		for (e = 0; e < s->di; e++)
-			inv[e] = 1.0 - in[e];
+			inv[e] = 1.0 - icx_powlike(in[e], s->ixpow);
 	}
 
 	/* If we have some sort of perceptual conversion */
@@ -273,8 +285,9 @@ pcpt_to_nLab(pcpt *s, double *out, double *in) {
 		else if (s->mlu != NULL) {
 			s->mlu->lookup(s->mlu, lab, inv);
 			icmXYZ2Lab(&icmD50, lab, lab);
-		} else 
+		} else { 
 			s->clu->dev_to_rLab(s->clu, lab, inv);
+		}
 
 #ifdef EMPH_NEUTRAL		/* Emphasise neutral axis, like CIE94 does */
 		{
@@ -283,13 +296,21 @@ pcpt_to_nLab(pcpt *s, double *out, double *in) {
 			c = sqrt(lab[1] * lab[1] + lab[2] * lab[2]);	/* Compute chromanance */
 
 //			c = 2.6624 / (1.0 + 0.013 * c);		/* Full strength scale factor */
-			c = 3.0 / (1.0 + 0.03 * c);		/* Full strength scale factor */
-			c = 1.0 + s->nemph * (c - 1.0);			/* Reduced strength scale factor */
+			c = 3.0 / (1.0 + 0.03 * c);			/* Full strength scale factor */
+			c = 1.0 + s->nemph * (c - 1.0);		/* Reduced strength scale factor */
 
 			lab[1] *= c;			/* scale a & b */
 			lab[2] *= c;
 		}
 #endif
+
+		/* Dark emphasis */
+		/* This doesn't actually match how demph is applied to device values... */
+		if (s->idemph < 1.0) {
+			double vv = lab[0];
+			lab[0] = 100.0 * pow(lab[0]/100.0, s->idemph);
+		}
+
 		/* Copy Lab values to output */
 		for (e = 0; e < (s->di < 3 ? s->di : 3); e++)
 			out[e] = lab[e];
@@ -611,7 +632,9 @@ inkmask xmask,			/* external xcolorants mask */
 inkmask nmask,			/* internal xcolorants mask */
 double *ilimit,			/* ink sum limit (scale 1.0) input and return, -1 if default */
 double *uilimit,		/* underlying ink sum limit (scale 1.0) input and return, -1 if default */
-double nemph			/* Neutral emphasis, 0.0 - 1.0. < 0.0 for default == CIE94 */
+double nemph,			/* Neutral emphasis, 0.0 - 1.0. < 0.0 for default == CIE94 */
+double demph,			/* Dark emphasis, 1.0 - 4.0. < 0.0 for default == none */
+double xpow				/* Extra device power, default = none */
 ) {
 	int e;
 	pcpt *s;
@@ -638,6 +661,14 @@ double nemph			/* Neutral emphasis, 0.0 - 1.0. < 0.0 for default == CIE94 */
 		nemph = NEMPH_DEFAULT;
 	s->nemph = nemph;
 
+	if (demph < 0.0)
+		demph = DEMPH_DEFAULT;
+	s->idemph = demph;
+
+	if (xpow < 0.0)
+		xpow = XPOW_DEFAULT;
+	s->ixpow = xpow;
+
 	/* See if we have a profile */
 	if (profName != NULL
 	 && profName[0] != '\000'
@@ -809,9 +840,11 @@ usage(int level, char *diag, ...) {
 	}
 	fprintf(stderr," -G               Generate good optimized points rather than Fast\n");
 	fprintf(stderr," -e patches       White test patches (default 4)\n");
+	fprintf(stderr," -B patches       Black test patches (default 4 Grey/RGB, else 0)\n");
 	fprintf(stderr," -s steps         Single channel steps (default grey 50, color 0)\n");
 	fprintf(stderr," -g steps         Grey axis RGB or CMY steps (default 0)\n");
 	fprintf(stderr," -m steps         Multidimensional device space cube steps (default 0)\n");
+	fprintf(stderr," -b steps         Multidimensional body centered cubic steps (default 0)\n");
 	fprintf(stderr," -f patches       Add iterative & adaptive full spread patches to total (default grey 0, color 836)\n");
 	fprintf(stderr,"                  Default is Optimised Farthest Point Sampling (OFPS)\n");
 	fprintf(stderr,"  -t              Use incremental far point for full spread\n");
@@ -830,7 +863,8 @@ usage(int level, char *diag, ...) {
 	fprintf(stderr," -p power         Optional power-like value applied to all device values.\n");
 	fprintf(stderr," -c profile       Optional device ICC or MPP pre-conditioning profile filename\n");
 	fprintf(stderr,"                  (Use \"none\" to turn off any conditioning)\n");
-	fprintf(stderr," -N emphasis      Degree of neutral axis patch concentration 0.0-1.0 (default %.2f)\n",NEMPH_DEFAULT);
+	fprintf(stderr," -N nemphasis     Degree of neutral axis patch concentration 0.0-1.0 (default %.2f)\n",NEMPH_DEFAULT);
+	fprintf(stderr," -V demphasis     Degree of dark region patch concentration 1.0-4.0 (default %.2f = none)\n",DEMPH_DEFAULT);
 	fprintf(stderr," -F L,a,b,rad     Filter out samples outside Lab sphere.\n");
 #ifdef VRML_DIAG
 	fprintf(stderr," -w               Dump diagnostic outfilel.wrl file (Lab locations)\n");
@@ -860,6 +894,8 @@ int dofilt(
 	return 0;
 }
 
+static double disprespt(cgats *pp, int p1, int p2);
+
 int main(int argc, char *argv[]) {
 	int i, j, k;
 	int fa, nfa, mfa;		/* current argument we're looking at */
@@ -873,10 +909,12 @@ int main(int argc, char *argv[]) {
 	char *ident;			/* Ink combination identifier (includes possible leading 'i') */
 	int good = 0;			/* 0 - fast, 1 = good */
 	int esteps = 4;			/* White color patches */
+	int Bsteps = -1;		/* Black color patches */
 	int ssteps = -1;		/* Single channel steps */
 	double xpow = 1.0;		/* Power to apply to all device values created */
 	int gsteps = 0;			/* Composite grey wedge steps */
 	int msteps = 0;			/* Regular grid multidimensional steps */
+	int bsteps = 0;			/* Regular body centered cubic grid multidimensional steps */
 	int fsteps = -1;		/* Fitted Multidimensional patches */
 	int uselat = 0;			/* Use incremental far point alg. for full spread points */
 	int userand = 0;		/* Use random for full spread points, 2 = perceptual */
@@ -890,6 +928,7 @@ int main(int argc, char *argv[]) {
 	double ilimit = -1.0;	/* Ink limit (scale 1.0) (default none) */
 	double uilimit = -1.0;	/* Underlying (pre-calibration, scale 1.0) ink limit */
 	double nemph = NEMPH_DEFAULT;
+	double demph = DEMPH_DEFAULT;
 	int filter = 0;			/* Filter values */
 	double filt[4] = { 50,0,0,0 };	
 	static char fname[MAXNAMEL+1] = { 0 };		/* Output file base name */
@@ -920,7 +959,7 @@ int main(int argc, char *argv[]) {
 	mfa = 1;        /* 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 */
+		if (argv[fa][0] == '-') {	/* Look for any flags */
 			char *na = NULL;		/* next argument after flag, null if none */
 
 			if (argv[fa][2] != '\000')
@@ -940,7 +979,7 @@ int main(int argc, char *argv[]) {
 				usage(0, "Usage requested");
 			}
 
-			else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') {
+			else if (argv[fa][1] == 'v') {
 				verb = 1;
 				if (na != NULL && na[0] >= '0' && na[0] <= '9') {
 					verb = atoi(na);
@@ -978,51 +1017,70 @@ int main(int argc, char *argv[]) {
 				good = 1;
 			}
 			/* White color patches */
-			else if (argv[fa][1] == 'e' || argv[fa][1] == 'E') {
+			else if (argv[fa][1] == 'e') {
 				int tt;
-				fa = nfa;
 				if (na == NULL) usage(0,"Expect argument after -e");
 				if ((tt = atoi(na)) >= 0)
 					esteps = tt;
+				fa = nfa;
 			}
-			/* Individual chanel steps */
-			else if (argv[fa][1] == 's' || argv[fa][1] == 'S') {
+			/* Black color patches */
+			else if (argv[fa][1] == 'B') {
 				int tt;
+				if (na == NULL) usage(0,"Expect argument after -B");
+				if ((tt = atoi(na)) >= 0)
+					Bsteps = tt;
 				fa = nfa;
+			}
+			/* Individual chanel steps */
+			else if (argv[fa][1] == 's') {
+				int tt;
 				if (na == NULL) usage(0,"Expect argument after -s");
 				if ((tt = atoi(na)) >= 0)
 					ssteps = tt;
+				fa = nfa;
 			}
 			/* RGB or CMY grey wedge steps */
 			else if (argv[fa][1] == 'g') {
 				int tt;
-				fa = nfa;
 				if (na == NULL) usage(0,"Expect argument after -g");
 				if ((tt = atoi(na)) >= 0)
 					gsteps = tt;
+				fa = nfa;
 			}
 			/* Multidimentional cube steps */
 			else if (argv[fa][1] == 'm') {
 				int tt;
-				fa = nfa;
 				if (na == NULL) usage(0,"Expect argument after -m");
 				if ((tt = atoi(na)) >= 0) {
 					msteps = tt;
 					if (msteps == 1)
 						msteps = 2;
 				}
+				fa = nfa;
+			}
+			/* Multidimentional body centered cube steps */
+			else if (argv[fa][1] == 'b') {
+				int tt;
+				if (na == NULL) usage(0,"Expect argument after -b");
+				if ((tt = atoi(na)) >= 0) {
+					bsteps = tt;
+					if (bsteps == 1)
+						bsteps = 2;
+				}
+				fa = nfa;
 			}
 			/* Full even spread Multidimentional patches */
 			else if (argv[fa][1] == 'f') {
 				int tt;
-				fa = nfa;
 				if (na == NULL) usage(0,"Expect argument after -f");
 				if ((tt = atoi(na)) >= 0)
 					fsteps = tt;
+				fa = nfa;
 			}
 
 			/* Use incremental far point algorithm for full spread */
-			else if (argv[fa][1] == 't' || argv[fa][1] == 'T') {
+			else if (argv[fa][1] == 't') {
 				uselat = 1;
 				userand = 0;
 				useqrand = 0;
@@ -1031,7 +1089,8 @@ int main(int argc, char *argv[]) {
 			}
 
 			/* Random requested */
-			else if (argv[fa][1] == 'r' || argv[fa][1] == 'R') {
+			else if (argv[fa][1] == 'r'
+				  || argv[fa][1] == 'R') {
 				uselat = 0;
 				if (argv[fa][1] == 'R')
 					userand = 2;
@@ -1043,7 +1102,8 @@ int main(int argc, char *argv[]) {
 			}
 
 			/* Space filling quasi-random requested */
-			else if (argv[fa][1] == 'q' || argv[fa][1] == 'Q') {
+			else if (argv[fa][1] == 'q'
+				  || argv[fa][1] == 'Q') {
 				uselat = 0;
 				userand = 0;
 				if (argv[fa][1] == 'Q')
@@ -1075,14 +1135,13 @@ int main(int argc, char *argv[]) {
 
 			/* Simplex grid angle */
 			else if (argv[fa][1] == 'a') {
-				fa = nfa;
 				if (na == NULL) usage(0,"Expect argument after -a");
 				simangle = atof(na);
+				fa = nfa;
 			}
 
 			/* Degree of iterative adaptation */
 			else if (argv[fa][1] == 'A') {
-				fa = nfa;
 				if (na == NULL) usage(0,"Expected argument to average deviation flag -A");
 				if (na[0] == 'p') {			/* (relative, for verification) */
 					perc_wght = atof(na+1);
@@ -1099,49 +1158,59 @@ int main(int argc, char *argv[]) {
 					if (dadapt < 0.0 || dadapt > 1.0)
 						usage(0,"Average Deviation argument %f must be between 0.0 and 1.0",dadapt);
 				}
+				fa = nfa;
 			}
 
 			/* Ink limit percentage */
-			else if (argv[fa][1] == 'l' || argv[fa][1] == 'L') {
+			else if (argv[fa][1] == 'l') {
 				double tt;
-				fa = nfa;
 				if (na == NULL) usage(0,"Expect argument after -l");
 				if ((tt = atof(na)) > 0.0)
 					uilimit = ilimit = 0.01 * tt;
+				fa = nfa;
 			}
 
 			/* Extra device power-like to use */
-			else if (argv[fa][1] == 'p' || argv[fa][1] == 'P') {
+			else if (argv[fa][1] == 'p') {
 				double tt;
-				fa = nfa;
 				if (na == NULL) usage(0,"Expect argument after -p");
 				if ((tt = atof(na)) > 0.0)
 					xpow = tt;
+				fa = nfa;
 			}
 
 			/* ICC profile for perceptual linearisation */
-			else if (argv[fa][1] == 'c' || argv[fa][1] == 'C') {
-				fa = nfa;
+			else if (argv[fa][1] == 'c') {
 				if (na == NULL) usage(0,"Expect argument after -c");
 				strncpy(pname,na,MAXNAMEL-1); pname[MAXNAMEL-1] = '\000';
+				fa = nfa;
 			}
 
 			/* Degree of neutral axis emphasis */
 			else if (argv[fa][1] == 'N') {
-				fa = nfa;
 				if (na == NULL) usage(0,"Expected argument to neutral emphasis flag -N");
 				nemph = atof(na);
 				if (nemph < 0.0 || nemph > 10.0)
 					usage(0,"Neautral weighting argument %f to '-N' is out of range",nemph);
+				fa = nfa;
+			}
+
+			/* Degree of dark region emphasis */
+			else if (argv[fa][1] == 'V') {
+				if (na == NULL) usage(0,"Expected argument to dark emphasis flag -V");
+				demph = atof(na);
+				if (demph < 1.0 || demph > 4.0)
+					usage(0,"Dark weighting argument %f to '-V' is out of range",demph);
+				fa = nfa;
 			}
 
 			/* Filter out samples outside given sphere */
 			else if (argv[fa][1] == 'F') {
-				fa = nfa;
 				if (na == NULL) usage(0,"Expect argument after -F");
 				if (sscanf(na, " %lf,%lf,%lf,%lf ",&filt[0], &filt[1], &filt[2], &filt[3]) != 4)
 					usage(0,"Argument to -F '%s' isn't correct",na);
 				filter = 1;
+				fa = nfa;
 			}
 
 #ifdef VRML_DIAG
@@ -1187,6 +1256,13 @@ int main(int argc, char *argv[]) {
 	stime = clock();
 
 	/* Implement some defaults */
+	if (Bsteps < 0) {
+		if (xmask == ICX_W || xmask == ICX_K || xmask == ICX_RGB || xmask == ICX_IRGB)
+			Bsteps = 4;
+		else
+			Bsteps = 0; 
+	}
+
 	if (di == 1) {
 		if (ssteps < 0)
 			ssteps = 50;
@@ -1201,22 +1277,22 @@ int main(int argc, char *argv[]) {
 
 	/* Do some sanity checking */
 	if (di == 1) {
-		if (ssteps == 0 && fsteps == 0 && msteps == 0)
+		if (ssteps == 0 && fsteps == 0 && msteps == 0 && bsteps == 0)
 			error ("Must have some Gray steps");
 		if (gsteps > 0) {
 			warning ("Composite grey steps ignored for monochrome output");
 			gsteps = 0;
 		}
 	} else if (di == 3) {
-		if (ssteps == 0 && fsteps == 0 && msteps == 0 && gsteps == 0)
+		if (ssteps == 0 && fsteps == 0 && msteps == 0 && bsteps == 0 && gsteps == 0)
 			error ("Must have some single or multi dimensional RGB or CMY steps");
 	} else {
-		if (ssteps == 0 && fsteps == 0 && msteps == 0 && gsteps == 0)
+		if (ssteps == 0 && fsteps == 0 && msteps == 0 && bsteps == 0 && gsteps == 0)
 			error ("Must have some single or multi dimensional steps");
 	}
 
 	/* Deal with ICC, MPP or fallback profile */
-	if ((pdata = new_pcpt(pname, xmask, nmask, &ilimit, &uilimit, nemph)) == NULL) {
+	if ((pdata = new_pcpt(pname, xmask, nmask, &ilimit, &uilimit, nemph, demph, xpow)) == NULL) {
 		error("Perceptual lookup object creation failed");
 	}
 
@@ -1231,6 +1307,10 @@ int main(int argc, char *argv[]) {
 	if (verb) {
 		printf("%s test chart\n",ident);
 
+		if (esteps > 0)
+			printf("White patches = %d\n",esteps);
+		if (Bsteps > 0)
+			printf("Black patches = %d\n",Bsteps);
 		if (ssteps > 0)
 			printf("Single channel steps = %d\n",ssteps);
 		if (gsteps > 0)
@@ -1239,6 +1319,8 @@ int main(int argc, char *argv[]) {
 			printf("Full spread patches = %d\n",fsteps);
 		if (msteps > 0)
 			printf("Multi-dimention cube steps = %d\n",msteps);
+		if (bsteps > 0)
+			printf("Multi-dimention body centered cube steps = %d\n",bsteps);
 		if (ilimit >= 0.0)
 			printf("Ink limit = %.1f%% (underlying %.1f%%)\n",ilimit * 100.0, uilimit * 100.0);
 		if (filter) {
@@ -1337,6 +1419,11 @@ int main(int argc, char *argv[]) {
 		pp->add_kword(pp, 0, "EXTRA_DEV_POW",buf, NULL);
 	}
 
+	if (demph > 1.0) {
+		sprintf(buf,"%f",demph);
+		pp->add_kword(pp, 0, "DARK_REGION_EMPHASIS",buf, NULL);
+	}
+
 	/* Only use optimsed full spread if <= 4 dimensions, else use ifarp */
 	if (di > 4 
 	 && userand == 0		/* Not other high D useful method */
@@ -1372,21 +1459,88 @@ int main(int argc, char *argv[]) {
 					val[e] = 0.0;			/* White is no colorant */
 				}
 			}
-	
+
 			/* Apply general filter */
 			if (filter && dofilt(pdata, filt, val))
 				continue;
 
 			sprintf(buf,"%d",id++);
 			ary[0].c = buf;
+
+			if (xmask == nmask) {
+				for (e = 0; e < di; e++)
+					ary[1 + e].d = 100.0 * val[e];
+			} else {
+				for (e = 0; e < di; e++)
+					ary[1 + e].d = 100.0 * (1.0 - val[e]);
+			}
+
 			pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
+			ary[1 + di + 0].d = 100.0 * XYZ[0];
+			ary[1 + di + 1].d = 100.0 * XYZ[1];
+			ary[1 + di + 2].d = 100.0 * XYZ[2];
+	
+			pp->add_setarr(pp, 0, ary);
+	
+			if (fxlist != NULL) {		/* Note in fixed list */
+				if (fxno >= fxlist_a) {
+					fxlist_a *= 2;
+					if ((fxlist = (fxpos *)realloc(fxlist, sizeof(fxpos) * fxlist_a)) == NULL)
+						error ("Failed to malloc fxlist");
+				}
+				for (e = 0; e < di; e++)
+					fxlist[fxno].p[e] = val[e];
+				fxlist[fxno].eloc = pp->t[0].nsets;
+				fxno++;
+			}
+		}
+	}
+
+	/* Black color patches */
+	if (Bsteps > 0)	{
+		int j, k, e;
+
+		for (j = k = 0; j < Bsteps; j++) {
+			double val[MXTD], XYZ[3];
+			cgats_set_elem ary[1 + MXTD + 3];
+	
+			if (nmask & ICX_ADDITIVE) {
+				for (e = 0; e < di; e++) {
+					val[e] = 0.0;			/* Black is no colorant */
+				}
+			} else {
+				for (e = 0; e < di; e++) {
+					val[e] = 1.0;			/* Black is full colorant */
+				}
+			}
+	
+			/* Apply general filter */
+			if (filter && dofilt(pdata, filt, val))
+				continue;
+
+			/* Do a simple ink limit */
+			if (uilimit < (double)di) {
+				double tot = 0.0;
+				for (e = 0; e < di; e++)
+					tot += val[e];
+				if (tot > uilimit) {
+					for (e = 0; e < di; e++)
+						val[e] *= uilimit/tot;
+				}
+			}
+
+			sprintf(buf,"%d",id++);
+			ary[0].c = buf;
+
 			if (xmask == nmask) {
 				for (e = 0; e < di; e++)
-					ary[1 + e].d = 100.0 * icx_powlike(val[e],xpow);
+					ary[1 + e].d = 100.0 * val[e];
 			} else {
 				for (e = 0; e < di; e++)
-					ary[1 + e].d = 100.0 * (1.0 - icx_powlike(val[e],xpow));
+					ary[1 + e].d = 100.0 * (1.0 - val[e]);
 			}
+
+			pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 			ary[1 + di + 0].d = 100.0 * XYZ[0];
 			ary[1 + di + 1].d = 100.0 * XYZ[1];
 			ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -1401,12 +1555,17 @@ int main(int argc, char *argv[]) {
 				}
 				for (e = 0; e < di; e++)
 					fxlist[fxno].p[e] = val[e];
+				fxlist[fxno].eloc = pp->t[0].nsets;
 				fxno++;
 			}
+			k++;
 		}
+		sprintf(buf,"%d",k);
+		pp->add_kword(pp, 0, "BLACK_COLOR_PATCHES",buf, NULL);
+	
 	}
 
-	/* Primary wedge steps */
+	/* Primary (single channel) wedge steps */
 	if (ssteps > 0)	{
 		sprintf(buf,"%d",ssteps);
 		pp->add_kword(pp, 0, "SINGLE_DIM_STEPS",buf, NULL);
@@ -1425,7 +1584,9 @@ int main(int argc, char *argv[]) {
 						val[e] = 0.0;
 				}
 
-				pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
+				/* Extra power and dark emphasis */
+				for (e = 0; e < di; e++)
+					val[e] = icx_powlike(val[e], xpow * demph);
 
 				/* See if it is already in the fixed list */
 				if (fxlist != NULL) {
@@ -1452,13 +1613,16 @@ int main(int argc, char *argv[]) {
 		
 					sprintf(buf,"%d",id++);
 					ary[0].c = buf;
+
 					if (xmask == nmask) {
 						for (e = 0; e < di; e++)
-							ary[1 + e].d = 100.0 * icx_powlike(val[e],xpow);
+							ary[1 + e].d = 100.0 * val[e];
 					} else {
 						for (e = 0; e < di; e++)
-							ary[1 + e].d = 100.0 * (1.0 - icx_powlike(val[e],xpow));
+							ary[1 + e].d = 100.0 * (1.0 - val[e]);
 					}
+
+					pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 					ary[1 + di + 0].d = 100.0 * XYZ[0];
 					ary[1 + di + 1].d = 100.0 * XYZ[1];
 					ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -1473,6 +1637,7 @@ int main(int argc, char *argv[]) {
 						}
 						for (e = 0; e < di; e++)
 							fxlist[fxno].p[e] = val[e];
+						fxlist[fxno].eloc = -1;
 						fxno++;
 					}
 				}
@@ -1514,16 +1679,17 @@ int main(int argc, char *argv[]) {
 					val[e] = 0.0;
 			}
 
+			/* Extra power and dark emphasis */
+			for (e = 0; e < di; e++)
+				val[e] = icx_powlike(val[e], xpow * demph);
+
 			/* Apply general filter */
 			if (filter && dofilt(pdata, filt, val))
 				addp = 0;
 
-			pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
-
-			/* Compute sum that includes affect of power */
+			/* Check if over ink limit */
 			for (sum = 0.0, e = 0; e < di; e++)
-				sum += icx_powlike(val[e], xpow);
-
+				sum += val[e];
 			if (sum > uilimit)
 				addp = 0;
 
@@ -1548,13 +1714,16 @@ int main(int argc, char *argv[]) {
 	
 				sprintf(buf,"%d",id++);
 				ary[0].c = buf;
+
 				if (xmask == nmask) {
 					for (e = 0; e < di; e++)
-						ary[1 + e].d = 100.0 * icx_powlike(val[e],xpow);
+						ary[1 + e].d = 100.0 * val[e];
 				} else {
 					for (e = 0; e < di; e++)
-						ary[1 + e].d = 100.0 * (1.0 - icx_powlike(val[e],xpow));
+						ary[1 + e].d = 100.0 * (1.0 - val[e]);
 				}
+
+				pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 				ary[1 + di + 0].d = 100.0 * XYZ[0];
 				ary[1 + di + 1].d = 100.0 * XYZ[1];
 				ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -1569,6 +1738,7 @@ int main(int argc, char *argv[]) {
 					}
 					for (e = 0; e < di; e++)
 						fxlist[fxno].p[e] = val[e];
+					fxlist[fxno].eloc = -1;
 					fxno++;
 				}
 			}
@@ -1594,16 +1764,17 @@ int main(int argc, char *argv[]) {
 			for (e = 0; e < di; e++)
 				val[e] = (double)gc[e]/(msteps-1);
 
+			/* Extra power and dark emphasis */
+			for (e = 0; e < di; e++)
+				val[e] = icx_powlike(val[e], xpow * demph);
+
 			/* Apply general filter */
 			if (filter && dofilt(pdata, filt, val))
 				addp = 0;
 
-			pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
-
-			/* Compute sum that includes affect of power */
+			/* Check if over ink limit */
 			for (sum = 0.0, e = 0; e < di; e++)
-				sum += icx_powlike(val[e], xpow);
-
+				sum += val[e];
 			if (sum > uilimit)
 				addp = 0;		/* Don't add patches over ink limit */
 
@@ -1630,13 +1801,16 @@ int main(int argc, char *argv[]) {
 
 				sprintf(buf,"%d",id++);
 				ary[0].c = buf;
+
 				if (xmask == nmask) {
 					for (e = 0; e < di; e++)
-						ary[1 + e].d = 100.0 * icx_powlike(val[e],xpow);
+						ary[1 + e].d = 100.0 * val[e];
 				} else {
 					for (e = 0; e < di; e++)
-						ary[1 + e].d = 100.0 * (1.0 - icx_powlike(val[e],xpow));
+						ary[1 + e].d = 100.0 * (1.0 - val[e]);
 				}
+
+				pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 				ary[1 + di + 0].d = 100.0 * XYZ[0];
 				ary[1 + di + 1].d = 100.0 * XYZ[1];
 				ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -1651,11 +1825,11 @@ int main(int argc, char *argv[]) {
 					}
 					for (e = 0; e < di; e++)
 						fxlist[fxno].p[e] = val[e];
+					fxlist[fxno].eloc = -1;
 					fxno++;
 				}
 			}
 
-	next_cpoint:;
 			/* Increment grid index and position */
 			for (j = 0; j < di; j++) {
 				gc[j]++;
@@ -1719,14 +1893,14 @@ int main(int argc, char *argv[]) {
 							if (addp) {
 								sprintf(buf,"%d",id++);
 								ary[0].c = buf;
-								pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 								if (xmask == nmask) {
 									for (e = 0; e < di; e++)
-										ary[1 + e].d = 100.0 * icx_powlike(val[e],xpow);
+										ary[1 + e].d = 100.0 * val[e];;
 								} else {
 									for (e = 0; e < di; e++)
-										ary[1 + e].d = 100.0 * (1.0 - icx_powlike(val[e],xpow));
+										ary[1 + e].d = 100.0 * (1.0 - val[e]);
 								}
+								pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 								ary[1 + di + 0].d = 100.0 * XYZ[0];
 								ary[1 + di + 1].d = 100.0 * XYZ[1];
 								ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -1741,6 +1915,7 @@ int main(int argc, char *argv[]) {
 									}
 									for (e = 0; e < di; e++)
 										fxlist[fxno].p[e] = val[e];
+									fxlist[fxno].eloc = -1;
 									fxno++;
 								}
 							}
@@ -1762,6 +1937,107 @@ int main(int argc, char *argv[]) {
 #endif /* ADDRECCLIPPOINTS */
 	}
 
+	/* Regular body centered cubic gridded Multi dimension steps */
+	if (bsteps > 0) {
+		int gc[MXTD];			/* Grid coordinate */
+		int pass = 0;			/* 0 = outer grid, 1 = inner grid */
+
+		sprintf(buf,"%d",bsteps);
+		pp->add_kword(pp, 0, "MULTI_DIM_BCC_STEPS",buf, NULL);
+
+		for (pass = 0; pass < 2; pass++) {
+
+			for (j = 0; j < di; j++)
+				gc[j] = 0;			/* init coords */
+				
+			for (;;) {	/* For all grid points */
+				double sum, val[MXTD], XYZ[3];
+				int addp, e;
+
+				addp = 1;			/* Default add the point */
+
+				for (e = 0; e < di; e++)
+					val[e] = (double)(pass * 0.5 + gc[e])/(bsteps-1);
+
+				/* Extra power and dark emphasis */
+				for (e = 0; e < di; e++)
+					val[e] = icx_powlike(val[e], xpow * demph);
+
+				/* Apply general filter */
+				if (filter && dofilt(pdata, filt, val))
+					addp = 0;
+
+				/* Check if over ink limit */
+				for (sum = 0.0, e = 0; e < di; e++)
+					sum += val[e];
+				if (sum > uilimit)
+					addp = 0;		/* Don't add patches over ink limit */
+
+				/* See if it is already in the fixed list */
+				if (addp && fxlist != NULL) { 
+					int k;
+					for (k = 0; k < fxno; k++) {
+						for (e = 0; e < di; e++) {
+							double tt;
+							tt = fabs(fxlist[k].p[e] - val[e]);
+							if (tt > MATCH_TOLL)
+								break;			/* Not identical */
+						}
+						if (e >= di)
+							break;				/* Was identical */
+					}
+					if (k < fxno)				/* Found an identical patch */
+						addp = 0;				/* Don't add the point */
+				}
+
+				/* Add patch to list if OK */
+				if (addp) {
+					cgats_set_elem ary[1 + MXTD + 3];
+
+					sprintf(buf,"%d",id++);
+					ary[0].c = buf;
+					if (xmask == nmask) {
+						for (e = 0; e < di; e++)
+							ary[1 + e].d = 100.0 * val[e];
+					} else {
+						for (e = 0; e < di; e++)
+							ary[1 + e].d = 100.0 * (1.0 - val[e]);
+					}
+
+					pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
+					ary[1 + di + 0].d = 100.0 * XYZ[0];
+					ary[1 + di + 1].d = 100.0 * XYZ[1];
+					ary[1 + di + 2].d = 100.0 * XYZ[2];
+
+					pp->add_setarr(pp, 0, ary);
+
+					if (fxlist != NULL) {		/* Note in fixed list */
+						if (fxno >= fxlist_a) {
+							fxlist_a *= 2;
+							if ((fxlist = (fxpos *)realloc(fxlist, sizeof(fxpos) * fxlist_a)) == NULL)
+								error ("Failed to malloc fxlist");
+						}
+						for (e = 0; e < di; e++)
+							fxlist[fxno].p[e] = val[e];
+						fxlist[fxno].eloc = -1;
+						fxno++;
+					}
+				}
+
+				/* Increment grid index and position */
+				for (j = 0; j < di; j++) {
+					gc[j]++;
+					if ((pass == 0 && gc[j] < bsteps)
+					 || (pass == 1 && gc[j] < (bsteps-1)))
+						break;	/* No carry */
+					gc[j] = 0;
+				}
+				if (j >= di)
+					break;		/* Done grid */
+			}
+		}
+	}
+
 	if (fsteps > fxno) { /* Top up with full spread (perceptually even) and other patch types */
 
 		/* Generate device random numbers. Don't check for duplicates */
@@ -1791,16 +2067,17 @@ int main(int argc, char *argv[]) {
 						val[e] = d_rand(0.0, 1.0);
 				}
 
+				/* Extra power and dark emphasis */
+				for (e = 0; e < di; e++)
+					val[e] = icx_powlike(val[e], xpow * demph);
+
 				/* Apply general filter */
 				if (filter && dofilt(pdata, filt, val))
 					continue;
 
-				pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
-
-				/* Compute sum that includes the affect of power */
+				/* Check if over ink limit */
 				for (sum = 0.0, e = 0; e < di; e++)
-					sum += icx_powlike(val[e], xpow);
-
+					sum += val[e];
 				if (sum > uilimit)
 					continue;
 
@@ -1808,11 +2085,13 @@ int main(int argc, char *argv[]) {
 				ary[0].c = buf;
 				if (xmask == nmask) {
 					for (e = 0; e < di; e++)
-						ary[1 + e].d = 100.0 * icx_powlike(val[e],xpow);
+						ary[1 + e].d = 100.0 * val[e];
 				} else {
 					for (e = 0; e < di; e++)
-						ary[1 + e].d = 100.0 * (1.0 - icx_powlike(val[e],xpow));
+						ary[1 + e].d = 100.0 * (1.0 - val[e]);
 				}
+
+				pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 				ary[1 + di + 0].d = 100.0 * XYZ[0];
 				ary[1 + di + 1].d = 100.0 * XYZ[1];
 				ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -1827,6 +2106,7 @@ int main(int argc, char *argv[]) {
 					}
 					for (e = 0; e < di; e++)
 						fxlist[fxno].p[e] = val[e];
+					fxlist[fxno].eloc = -1;
 					fxno++;
 				}
 
@@ -1917,16 +2197,14 @@ int main(int argc, char *argv[]) {
 				if (filter && dofilt(pdata, filt, val))
 					continue;
 
-				pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
-
-				/* Do a simple ink limit that include the effect of xpow */
+				/* Do a simple ink limit */
 				if (uilimit < (double)di) {
 					double tot = 0.0;
 					for (e = 0; e < di; e++)
-						tot += icx_powlike(val[e],xpow);
+						tot += val[e];
 					if (tot > uilimit) {
 						for (e = 0; e < di; e++)
-							val[e] = icx_powlike(icx_powlike(val[e],xpow) * uilimit/tot, 1.0/xpow);
+							val[e] = val[e] * uilimit/tot;
 					}
 				}
 
@@ -1934,11 +2212,13 @@ int main(int argc, char *argv[]) {
 				ary[0].c = buf;
 				if (xmask == nmask) {
 					for (e = 0; e < di; e++)
-						ary[1 + e].d = 100.0 * icx_powlike(val[e],xpow);
+						ary[1 + e].d = 100.0 * val[e];
 				} else {
 					for (e = 0; e < di; e++)
-						ary[1 + e].d = 100.0 * (1.0 - icx_powlike(val[e],xpow));
+						ary[1 + e].d = 100.0 * (1.0 - val[e]);
 				}
+
+				pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 				ary[1 + di + 0].d = 100.0 * XYZ[0];
 				ary[1 + di + 1].d = 100.0 * XYZ[1];
 				ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -1953,6 +2233,7 @@ int main(int argc, char *argv[]) {
 					}
 					for (e = 0; e < di; e++)
 						fxlist[fxno].p[e] = val[e];
+					fxlist[fxno].eloc = -1;
 					fxno++;
 				}
 			}
@@ -1960,6 +2241,255 @@ int main(int argc, char *argv[]) {
 		}
 	}
 
+	/* Even the location of marked patches into sequence */
+	{
+		int ii, p1, p2, t1;
+
+		/* For each patch to be dispersed */
+		for (ii = 0; ii < fxno; ii++) {
+			if (fxlist[ii].eloc >= 0) {
+				p1 = fxlist[ii].eloc;
+
+				for (k = 0; k < 10; k++) {		/* Retry 10 times */
+
+					/* Pick a random patch to exchange it with */
+					p2 = i_rand(0, pp->t[0].nsets-1);
+	
+					/* Check it isn't one of our patches to be dispersed */
+					for (i = 0; i < fxno; i++) {
+						if (fxlist[i].eloc == p2)
+							break;
+					}
+					if (i < fxno)
+						continue;		/* Try another patch to exchange with */
+
+					/* Swap */
+					for (j = 1; j < (1 + di + 3); j++) {
+						double tt = *((double *)pp->t[0].fdata[p1][j]);
+						*((double *)pp->t[0].fdata[p1][j]) = *((double *)pp->t[0].fdata[p2][j]);
+						*((double *)pp->t[0].fdata[p2][j]) = tt;
+					}
+					fxlist[ii].eloc = p2;
+
+					break;
+				}
+			}
+		}
+	}
+
+	/* If this seems to be for a CRT, optimise the patch order to minimise the */
+	/* response time delays */
+	if (nmask == ICX_RGB && pp->t[0].nsets > 1) {
+		int npat = pp->t[0].nsets;
+		char *nm;						/* Don't move array */
+		double udelay, *delays, adelay;
+		double temp, trate;	/* Annealing temperature & rate */
+		double tstart, tend;/* Annealing chedule range */
+
+		if ((nm = (char *)malloc(sizeof(char) * npat)) == NULL)
+			error ("Failed to malloc nm array");
+		if ((delays = (double *)malloc(sizeof(double) * npat)) == NULL)
+			error ("Failed to malloc delay array");
+
+		/* Set nm[] to mark patches that shouldn't be moved */
+		for (i = 0; i < npat; i++)
+			nm[i] = 0;
+		for (i = 0; i < fxno; i++) {
+			if (fxlist[i].eloc >= 0)
+				nm[fxlist[i].eloc] = 1;
+		}
+
+#ifdef NEVER
+		/* Randomly shuffle patches */
+		{
+			int p1, p2;
+
+			for (p1 = 0; p1 < npat; p1++) {
+	
+				if (nm[p1])
+					continue;
+
+				p2 = i_rand(0, npat-1);
+				if (nm[p2])
+					continue;
+				for (j = 1; j < (1 + di + 3); j++) {
+					double tt = *((double *)pp->t[0].fdata[p1][j]);
+					*((double *)pp->t[0].fdata[p1][j]) = *((double *)pp->t[0].fdata[p2][j]);
+					*((double *)pp->t[0].fdata[p2][j]) = tt;
+				}
+			}
+		}
+#endif
+#ifdef NEVER
+		/* Simple sort by brightness */
+		{
+			int p1, p2;
+			double rgb1, rgb2;
+
+			for (p1 = 0; p1 < (npat-1); p1++) {
+	
+				if (nm[p1])
+					continue;
+
+				rgb1 = pow(*((double *)pp->t[0].fdata[p1][1 + 0]), 2.2)
+				     + pow(*((double *)pp->t[0].fdata[p1][1 + 1]), 2.2)
+				     + pow(*((double *)pp->t[0].fdata[p1][1 + 2]), 2.2);
+	
+				for (p2 = p1 + 1; p2 < npat; p2++) {
+
+					if (nm[p2])
+						continue;
+
+					rgb2 = pow(*((double *)pp->t[0].fdata[p2][1 + 0]), 2.2)
+					     + pow(*((double *)pp->t[0].fdata[p2][1 + 1]), 2.2)
+					     + pow(*((double *)pp->t[0].fdata[p2][1 + 2]), 2.2);
+		
+					if (rgb2 < rgb1) {
+						for (j = 1; j < (1 + di + 3); j++) {
+							double tt = *((double *)pp->t[0].fdata[p1][j]);
+							*((double *)pp->t[0].fdata[p1][j]) = *((double *)pp->t[0].fdata[p2][j]);
+							*((double *)pp->t[0].fdata[p2][j]) = tt;
+						}
+						rgb1 = rgb2;
+					}
+				}
+			}
+		}
+#endif
+		/* Compute the current overall update delay */
+		udelay = 0.0;
+		for (i = 1; i < npat; i++) {
+			double xdelay;
+
+			xdelay = disprespt(pp, i-1, i);
+
+			delays[i] = xdelay;
+//printf("~1 delay[%d] = %f\n",i,xdelay);
+			udelay += xdelay;
+		}
+
+		if (verb)
+			printf("Extra display response delay = %f sec., optimizing....\n",udelay);
+
+		{
+			int nchunks, chsize;
+			int chstart, chend;
+
+			if (verb)
+				printf("%c%2d%%",cr_char,0); fflush(stdout);
+
+			/* We'll do this in chunks of 500 to make it linear time overall, */
+			/* at the cost of the best possible optimisation. */
+			nchunks = (int)ceil(npat/500.0); 
+			chsize = (int)ceil(npat/nchunks);
+			for (chstart = 0; chstart < npat; chstart += chsize) {
+				int p1, p2, bp2;
+				double p1d, p2d, p1d1, p2d1;
+				double p1nd, p2nd, p1nd1, p2nd1;
+				double tdelay, de;
+				int noswapped;
+
+				chend = chstart + chsize+2;
+				if (chend > npat)
+					chend = npat;
+				noswapped = chend - chstart;
+//printf("~1 chstart %d, chend %d, size %d\n",chstart,chend, chend - chstart);
+
+				/* While we are still improving, and the improvement was significant */
+				for (;noswapped > 5;) {
+					noswapped = 0;
+
+					for (p1 = chstart + 1; p1 < chend; p1++) {
+						if (nm[p1])
+							continue;
+
+						p1d = delays[p1];
+
+						/* Locate the patch ahead of us that is best to swap with */
+						bp2 = -1;
+						for (p2 = p1 + 2; p2 < chend; p2++) {
+
+							if (nm[p2])
+								continue;
+
+							/* Compute effect of a swap on the total delay */
+							p2d = delays[p2];
+							p1nd = disprespt(pp, p2-1, p1);
+							p2nd = disprespt(pp, p1-1, p2);
+							p1d1 = p1nd1 = 0.0;
+							if ((p1+1) < chend) {
+								p1d1 = delays[p1+1];
+								p1nd1 = disprespt(pp, p2, p1+1);
+							}
+							p2d1 = p2nd1 = 0.0;
+							if ((p2+1) < chend) {
+								p2d1 = delays[p2+1];
+								p2nd1 = disprespt(pp, p1, p2+1);
+							}
+			
+							tdelay = udelay - p1d - p2d - p1d1 - p2d1 + p1nd + p2nd + p1nd1 + p2nd1;
+
+							if (tdelay < udelay) {
+								bp2 = p2;
+							}
+						}
+						if (bp2 < 0) {
+							continue;
+						} 
+
+						noswapped++;
+
+						p2 = bp2;
+
+						p2d = delays[p2];
+						p1nd = disprespt(pp, p2-1, p1);
+						p2nd = disprespt(pp, p1-1, p2);
+						p1d1 = p1nd1 = 0.0;
+						if ((p1+1) < chend) {
+							p1d1 = delays[p1+1];
+							p1nd1 = disprespt(pp, p2, p1+1);
+						}
+						p2d1 = p2nd1 = 0.0;
+						if ((p2+1) < chend) {
+							p2d1 = delays[p2+1];
+							p2nd1 = disprespt(pp, p1, p2+1);
+						}
+
+						tdelay = udelay - p1d - p2d - p1d1 - p2d1 + p1nd + p2nd + p1nd1 + p2nd1;
+
+						/* Swap the values */
+						udelay = tdelay;
+						delays[p2] = p1nd;
+						delays[p1] = p2nd;
+						if (p1 < (chend-1))
+							delays[p1+1] = p1nd1;
+						if (p2 < (chend-1))
+							delays[p2+1] = p2nd1;
+
+						for (j = 1; j < (1 + di + 3); j++) {
+							double tt = *((double *)pp->t[0].fdata[p1][j]);
+							*((double *)pp->t[0].fdata[p1][j]) = *((double *)pp->t[0].fdata[p2][j]);
+							*((double *)pp->t[0].fdata[p2][j]) = tt;
+						}
+//printf("~1 swaping %d and %d, udelay %f\n",p1,p2,udelay);
+					}
+//printf("~1 udelay %f\n",udelay);
+					if (verb) {
+						printf("%c%2d%%",cr_char,(int)(100.0 * (chend-1 - noswapped)/(npat-1.0)));
+						fflush(stdout);
+					}
+				}
+			}
+			if (verb)
+				printf("%c%2d%%",cr_char,100); fflush(stdout);
+		}
+		if (verb)
+			printf("\nOptimised display response delay = %f sec.\n",udelay);
+
+		free(delays);
+		free(nm);
+	}
+
 	/* Use ofps to measure the stats of the points */
 	/* Note that if new_ofps() fails it will exit() */
 	if (verb > 1
@@ -2001,11 +2531,6 @@ int main(int argc, char *argv[]) {
 
 			/* Lookup device values for target density */
 			pdata->den_to_dev(pdata, val, den);	
-			pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
-
-			/* Apply extra power */
-			for (e = 0; e < di; e++)
-				val[e] = icx_powlike(val[e], xpow);
 
 			/* Do a simple ink limit */
 			if (uilimit < (double)di) {
@@ -2027,6 +2552,7 @@ int main(int argc, char *argv[]) {
 				for (e = 0; e < di; e++)
 					ary[1 + e].d = 100.0 * (1.0 - val[e]);
 			}
+			pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 			ary[1 + di + 0].d = 100.0 * XYZ[0];
 			ary[1 + di + 1].d = 100.0 * XYZ[1];
 			ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -2090,18 +2616,16 @@ int main(int argc, char *argv[]) {
 					val[0] = val[1] = val[2] = 0.5;
 			}
 			
+			/* Apply extra power to device values (??) */
+			for (e = 0; e < di; e++)
+				val[e] = icx_powlike(val[e], xpow);
+
 			/* If target space isn't something we recognise, convert it */
 			if (ftarg != NULL) {
 				ftarg->dev_to_rLab(ftarg, lab, val);
 				pdata->rLab_to_dev(pdata, val, lab);	
 			}
 
-			pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
-
-			/* Apply extra power */
-			for (e = 0; e < di; e++)
-				val[e] = icx_powlike(val[e], xpow);
-
 			/* Do a simple ink limit */
 			if (uilimit < (double)di) {
 				double tot = 0.0;
@@ -2121,6 +2645,8 @@ int main(int argc, char *argv[]) {
 				for (e = 0; e < di; e++)
 					ary[1 + e].d = 100.0 * (1.0 - val[e]);
 			}
+
+			pdata->dev_to_XYZ(pdata, XYZ, val);		/* Add expected XYZ */
 			ary[1 + di + 0].d = 100.0 * XYZ[0];
 			ary[1 + di + 1].d = 100.0 * XYZ[1];
 			ary[1 + di + 2].d = 100.0 * XYZ[2];
@@ -2183,6 +2709,7 @@ int main(int argc, char *argv[]) {
 		rad = 15.0/pow(nsets, 1.0/(double)(di <= 3 ? di : 3));
 
 		for (i = 0; i < nsets; i++) {
+
 			/* Re-do any inversion before using dev_to_rLab() */
 			if (xmask == nmask) {
 				for (j = 0; j < di; j++)
@@ -2193,10 +2720,11 @@ int main(int argc, char *argv[]) {
 					idev[j] = 1.0 - dev[j];
 				}
 			}
+
 			pdata->dev_to_rLab(pdata, Lab, idev);
 			wrl->Lab2RGB(wrl, col, Lab);
 
-			/* Fudge device locations into Lab space */
+			/* Fudge device locations into "Lab" space */
 			Lab[0] = 100.0 * dev[0];
 			Lab[1] = 100.0 * dev[1] - 50.0;
 			Lab[2] = 100.0 * dev[2] - 50.0;
@@ -2217,6 +2745,47 @@ int main(int argc, char *argv[]) {
 	return 0;
 }
 
+/* - - - - - - - - - - - - - - - - - - - - - - - - - -  */
+
+/* Compte the display response time */
+static double disprespt(cgats *pp, int p1, int p2) {
+	double kr, kf;
+	double orgb[3], rgb[3];
+	double xdelay = 0.0;
+	int j;
+
+	kr = DISPLAY_RISE_TIME/log(1 - 0.9);	/* Exponent constant */
+	kf = DISPLAY_FALL_TIME/log(1 - 0.9);	/* Exponent constant */
+
+	orgb[0] = *((double *)pp->t[0].fdata[p1][1 + 0]) / 100.0;
+	orgb[1] = *((double *)pp->t[0].fdata[p1][1 + 1]) / 100.0;
+	orgb[2] = *((double *)pp->t[0].fdata[p1][1 + 2]) / 100.0;
+
+	rgb[0] = *((double *)pp->t[0].fdata[p2][1 + 0]) / 100.0;
+	rgb[1] = *((double *)pp->t[0].fdata[p2][1 + 1]) / 100.0;
+	rgb[2] = *((double *)pp->t[0].fdata[p2][1 + 2]) / 100.0;
+
+	for (j = 0; j < 3; j++) {
+		double el, dl, n, t;
+
+		el = pow(rgb[j], 2.2);
+		dl = el - pow(orgb[j], 2.2);	/* Change in level */
+		if (fabs(dl) > 0.01) {		/* More than 1% change in level */
+			n = DISPLAY_SETTLE_AIM * el;
+			if (n < DISPLAY_ABS_AIM)
+				n = DISPLAY_ABS_AIM;
+			if (dl > 0.0)
+				t = kr * log(n/dl);
+			else
+				t = kf * log(n/-dl);
+
+			if (t > xdelay)
+				xdelay = t;
+		}
+	}
+	return xdelay;
+}
+