Imported Upstream version 1.6.3upstream/1.6.3

4 files changed, 866 insertions, 130 deletions

diff --git a/target/ofps.c b/target/ofps.c
index fdea8cb..e485035 100644
--- a/target/ofps.c
+++ b/target/ofps.c
@@ -32,7 +32,7 @@
 
 	Some profiles are too rough, and slow/stall vertex placement.
 	Reducing the cache grid and/or smoothing the rspl values
-	ay mitigate this to some degree, and make this more robust ??
+	may mitigate this to some degree, and make this more robust ??
 
  */
 
@@ -45,7 +45,7 @@
 
 	We then initially add sampling points at the largest estimated error
 	verticies of the veronoi natural neighbourhood.
-	This gives us an optimal distribution measuring in mestimated position
+	This gives us an optimal distribution measuring in estimated position
 	error within a tollerance of 2:1
 	
 	We then iteratively improve the distribution of point nodes by
@@ -84,7 +84,7 @@
 	all possible device channel values. This also applies at
 	higher dimensions (ie. the CMYK values exploring response
 	to different K values doesn't spread the CMY values
-	evenly appart.)
+	evenly apart.)
 
 	Stratification seems to be somewhat at odds with the primary goal
 	of minimizing the maximum estimated error from any point in the
diff --git a/target/printtarg.c b/target/printtarg.c
index 25d889a..3ad3023 100644
--- a/target/printtarg.c
+++ b/target/printtarg.c
@@ -17,12 +17,20 @@
 
 	TTBD:
 
+	Add independent w & h patch scaling option.
+
+	Allow scaling minimum leading/trailing white space.
+
+	Add option to omit labelling.
+
 	Add "single pixel patch" mode, for pure digital processing for
 	abstract profile creation.
 
 	Add -h2 flag for Munki for super high-res chart ?
 	Note:   i1Pro:  Illum spot: 3.5mm Aperture: 4.5mm, Physical aperture: 4.55mm
+
             Munki:  Illum spot: 8.0mm Aperture: 6.0mm, Physical aperture: 7.63mm
+			Min patch size is 6mm x 6mm, below that increases delta E.
 
 	Add an option that allows including a scale gauge, to detect
 	accidental re-scaling.
@@ -31,6 +39,8 @@
 	rather than at the trailing edges.
 
 	Add direct PDF support, including NChannel output.
+
+	Add option to apply a scale to counteract the Adobe utility problem.
 */
 
 /* This program generates a PostScript or TIFF print target file, */
@@ -98,7 +108,7 @@
  *    Improve EPS support to add a preview to each eps file.
  */
 
-#undef DEBUG
+#undef DEBUG			/* Print edge details to stderr */
 #undef FORCEN			/* For testing, force DeviceN */
 #define DEN_COMPRESS	/* Compress density estimates > 1.0 */
 						/* - this biases it towards white spacers */
@@ -526,6 +536,7 @@ trend *new_ps_trend(
 	int nmask,				/* Non zero if we are doing a DeviceN chart */
 	double pw, double ph,	/* Page width and height in mm */
 	int eps,				/* EPS flag */
+	int nocups,				/* NZ to supress cups job ticket */
 	int rand,				/* randomize */
 	int rstart				/* Random start number/chart ID */
 ) {
@@ -595,6 +606,8 @@ trend *new_ps_trend(
 		fprintf(s->of,"%%%%PageOrder: Ascend\n");
 		fprintf(s->of,"%%%%BoundingBox: %d %d %d %d\n",0,0,ipw-1,iph-1);
 		fprintf(s->of,"%%%%Orientation: Portrait\n");		/* Rows are always virtical */
+		if (!nocups)
+			fprintf(s->of,"%%cupsJobTicket: cups-disable-cmm\n");
 		fprintf(s->of,"%%%%EndComments\n");
 		fprintf(s->of,"\n");
 		if (!eps) {
@@ -1733,6 +1746,7 @@ int hflag,			/* Spectroscan/Munki high density modified */
 int verb,			/* Verbose flag */
 int scanc,			/* Scan compatible bits, 1 = .cht gen, 2 = wide first row */
 int oft,			/* PS/EPS/TIFF select (0,1,2) */
+int nocups,			/* NZ to supress cups job ticket in PS/EPS */
 depth2d tiffdpth,	/* TIFF pixel depth */
 double tiffres,		/* TIFF resolution in DPI */
 int ncha,			/* flag, use nchannel alpha */
@@ -2438,7 +2452,7 @@ int *p_npat			/* Return number of patches including padding */
 					if (oft == 0) {	/* PS */
 						if (flags & IS_FPIF) {					/* First page */
 							sprintf(psname,"%s.ps",bname);
-							if ((tro = new_ps_trend(psname,npages,nmask,pw,ph,oft,rand,rstart)) == NULL)
+							if ((tro = new_ps_trend(psname,npages,nmask,pw,ph,oft,nocups,rand,rstart)) == NULL)
 								error ("Unable to create output rendering object file '%s'",psname);
 							if (verb)
 								printf("Creating file '%s'\n",psname);
@@ -2448,7 +2462,7 @@ int *p_npat			/* Return number of patches including padding */
 							sprintf(psname,"%s_%02d.eps",bname,pif);
 						else
 							sprintf(psname,"%s.eps",bname);
-						if ((tro = new_ps_trend(psname,npages,nmask,pw,ph,oft,rand,rstart)) == NULL)
+						if ((tro = new_ps_trend(psname,npages,nmask,pw,ph,oft,rand,rstart,nocups)) == NULL)
 							error ("Unable to create output rendering object file '%s'",psname);
 						if (verb)
 							printf("Creating file '%s'\n",psname);
@@ -2910,6 +2924,7 @@ void usage(char *diag, ...) {
 	fprintf(stderr," -M margin       Set a page margin in mm and include it in TIFF\n");       
 	fprintf(stderr," -P              Don't limit strip length\n");       
 	fprintf(stderr," -L              Suppress any left paper clip border\n");       
+	fprintf(stderr," -U              Suppress CUPS cupsJobTicket: cups-disable-cmm in PS & EPS files\n");
 	fprintf(stderr," -p size         Select page size from:\n");
 	for (pp = psizes; pp->name != NULL; pp++)
 		fprintf(stderr,"                 %-8s [%.1f x %.1f mm]%s\n", pp->name, pp->w, pp->h,
@@ -2932,6 +2947,7 @@ char *argv[];
 	int rand = 1;
 	int qbits = 0;			/* Quantization bits */
 	int oft = 0;			/* Ouput File type, 0 = PS, 1 = EPS , 2 = TIFF */
+	int nocups = 0;			/* Supress CUPS PS/EPS job ticket */
 	depth2d tiffdpth = bpc8_2d;	/* TIFF pixel depth */
 	double tiffres = 100.0;	/* TIFF resolution in DPI */
 	int ncha = 0;			/* flag, use nchannel alpha */
@@ -2996,7 +3012,8 @@ char *argv[];
 		usage("Not enough arguments");
 
 #ifdef DEBUG
-	printf("target: DEBUG is #defined\n");
+# pragma message("######### printtarg DEBUG is #define ########")
+	fprintf(stderr,"target: DEBUG is #defined\n");
 #endif
 
 	/* Find the default paper size */
@@ -3240,6 +3257,11 @@ char *argv[];
 				nolpcbord = 1;
 			}
 
+			/* Suppress CUPS job ticket */
+			else if (argv[fa][1] == 'U') {
+				nocups = 1;
+			}
+
 			/* Page size */
 			else if (argv[fa][1] == 'p') {
 				fa = nfa;
@@ -3646,12 +3668,12 @@ char *argv[];
 		pcol = pcold;		/* Density spacer alues */
 
 	
-	sprintf(label, "Argyll Color Management System - Test chart \"%s\" (%s %d) %s",
+	sprintf(label, "ArgyllCMS - Chart \"%s\" (%s %d) %s",
 	               psname, rand ? "Random Start" : "Chart ID", rstart, atm);
 	generate_file(itype, psname, cols, npat, applycal ? cal : NULL, label,
 	            pap != NULL ? pap->w : cwidth, pap != NULL ? pap->h : cheight,
 	            marg, nosubmarg, nollimit, nolpcbord, rand, rstart, saix, paix,	ixord,
-	            pscale, sscale, hflag, verb, scanc, oft, tiffdpth, tiffres, ncha, tiffdith,
+	            pscale, sscale, hflag, verb, scanc, oft, nocups, tiffdpth, tiffres, ncha, tiffdith,
 	            tiffcomp, spacer, nmask, altrep, pcol, wp,
 	            &sip, &pis, &plen, &glen, &tlen, &nppat);
 
@@ -3750,6 +3772,7 @@ char *argv[];
 /* A half edge structure */
 /* coordinate origin is top left */
 struct _hedge {
+	int ix;			/* Index for debug id */
 	double rgb[3];	/* Color this half edge transitions to */
 	int negh;		/* 1 if this is a -ve major coordinate side half edge, 0 otherwise */
 	double mj;		/* Major coordinate offset (ie. X coord for vertical edge) */
@@ -3761,7 +3784,7 @@ struct _hedge {
 /* A patch identifier */
 /* coordinate origin is top left */
 struct _patch {
-	char id[20];	/* ID string, Zeri length if a diagnostic rectangle */
+	char id[20];	/* ID string, Zero length if a diagnostic rectangle */
 	double xo;		/* Location of the rectangle origin (bottom left ???) */
 	double yo;
 	double w;		/* Size of the patch */
@@ -3799,7 +3822,7 @@ struct {
 
 	/* Raw half edge lists, [vertical, horizontal] */
 	int nhe[2];
-	hedge *he[2];
+	hedge *he[2];	/* Pointer to start of hedhe linked list */
 
 	/* Patch identity information */
 	int npatches;
@@ -3942,33 +3965,173 @@ double y
 /* .cht file. */
 void et_write(char *fname, col *cols, int *rix, int si, int ei) {
 	FILE *of;
-	hedge *ep;
+	hedge *ep0, *ep1, *epe;
 	int i, h;
 
 //printf("~1 et has %d vertical and %d horizontal half edges\n", et.nhe[0], et.nhe[1]);
 //printf("~1 et has %d patches\n", et.npatches);
 
+	/* Do X then Y */
 	for (h = 0; h < 2; h++) {
+
 		/* Create sorted list of vertical half edges */
 		if ((et.she[h] = (hedge **)malloc(sizeof(patch*) * et.nhe[h])) == NULL)
 			error("Malloc of array of vertical halfedge pointers failed");
 	
-		for (ep = et.he[h], i = 0; i < et.nhe[h]; i++, ep = ep->next)
-			et.she[h][i] = ep;
+		for (ep0 = et.he[h], i = 0; i < et.nhe[h]; i++, ep0 = ep0->next)
+			et.she[h][i] = ep0;
 	
 		/* Sort helf edges by their X location, then their Y0 location */
 #define HEAP_COMPARE(A,B) (fabs(A->mj - B->mj) < 1e-6 ? A->mi0 < B->mi0 : A->mj < B->mj)
 		HEAPSORT(hedge *, et.she[h], et.nhe[h]);
 #undef HEAP_COMPARE
 
-#ifdef NEVER
-for (i = 0; i < et.nhe[h]; i++) {
-printf("%s %d at %c = %f from %c = %f to %f\n",
-h == 0 ? "Vert" : "Horiz", i,
-h == 0 ? 'X' : 'Y', et.she[h][i]->mj,
-h == 0 ? 'Y' : 'X', et.she[h][i]->mi0, et.she[h][i]->mi1);
-}
-#endif /* NEVER */
+		/* Re-create the linked list in sorted order */
+		et.he[h] = NULL;
+		for (i = et.nhe[h]-1; i >= 0; i--) {
+			et.she[h][i]->next = et.he[h];
+			et.he[h] = et.she[h][i];
+			et.he[h]->ix = i;
+		}
+		
+		free(et.she[h]);
+		et.she[h] = NULL;
+
+#ifdef DEBUG
+		fprintf(stderr,"Sorted %s half edges:\n",h ? "Vertical" : "Horizontal");
+		for (ep0 = et.he[h]; ep0 != NULL; ep0 = ep0->next) {
+			fprintf(stderr,"%s %d at %c = %f from %c = %f to %f\n",
+			h == 0 ? "Vert" : "Horiz", ep0->ix,
+			h == 0 ? 'X' : 'Y', ep0->mj,
+			h == 0 ? 'Y' : 'X', ep0->mi0, ep0->mi1);
+		}
+#endif /* DEBUG */
+
+		/* Do a first pass to locate and split any part overlapping pairs of edges */
+		for (ep0 = et.he[h]; ep0 != NULL; ep0 = epe) {
+
+			/* Locate the end of the half edges at the same position */
+			for (epe = ep0; epe != NULL; epe = epe->next) {
+				if (epe == NULL || fabs(ep0->mj - epe->mj) > 1e-6)
+					break;
+			}
+
+#ifdef DEBUG
+			fprintf(stderr,"Doing group from %d to %d\n",ep0->ix, epe ? epe->ix : -1);
+#endif
+
+			/* Look for overlapping half edges, and split them up so */
+			/* there are no overlaps */
+			for (; ep0 != epe; ep0 = ep0->next) {
+
+				for (ep1 = ep0->next; ep1 != epe; ep1 = ep1->next) {
+
+					if (ep1->mi0 > ep0->mi1)	/* Out of range for overlap */
+						break;
+
+					/* If partial overlap */
+					if (ep0->mi0 < (ep1->mi0-1e-6)
+					 && ep0->mi1 > (ep1->mi0+1e-6)
+					 && ep0->mi1 < (ep1->mi1-1e-6)) {
+						hedge *ep0b, *ep1b;
+						
+#ifdef DEBUG
+						fprintf(stderr,"Half edges partial overlap:\n");
+						fprintf(stderr,"i = %d, j = %d\n",ep0->ix,ep1->ix);
+						fprintf(stderr,"%s %d at %c = %f from %c = %f to %f, half %s\n",
+						h == 0 ? "Vert" : "Horiz", ep0->ix,
+						h == 0 ? 'X' : 'Y', ep0->mj,
+						h == 0 ? 'Y' : 'X', ep0->mi0, ep0->mi1,
+						ep0->negh ? "Neg" : "Pos");
+						fprintf(stderr,"%s %d at %c = %f from %c = %f to %f, half %s\n",
+						h == 0 ? "Vert" : "Horiz", ep1->ix,
+						h == 0 ? 'X' : 'Y', ep1->mj,
+						h == 0 ? 'Y' : 'X', ep1->mi0, ep1->mi1,
+						ep1->negh ? "Neg" : "Pos");
+#endif
+						/* Split up the two edges so that we have four edges */
+						
+						if ((ep0b = (hedge *)calloc(sizeof(hedge), 1)) == NULL)
+							error("Malloc of half edge structure failed");
+						memcpy(ep0b, ep0, sizeof(hedge));
+
+						if ((ep1b = (hedge *)calloc(sizeof(hedge), 1)) == NULL)
+							error("Malloc of half edge structure failed");
+						memcpy(ep1b, ep1, sizeof(hedge));
+
+						ep0b->mi0 = ep1->mi0;
+						ep1b->mi1 = ep0->mi1;
+						ep0->mi1 = ep1->mi0;
+						ep1->mi0 = ep0->mi1;
+
+						/* Insert them in order into linked list */
+						ep1b->next = ep1;
+						ep0b->next = ep1b;
+						ep0->next = ep0b;
+
+						et.nhe[h] += 2;
+
+					/* If full overlap */
+					} else if (ep0->mi0 < (ep1->mi0-1e-6)
+					 && ep0->mi1 > (ep1->mi1+1e-6)) {
+						hedge *ep0b, *ep0c;
+						
+#ifdef DEBUG
+						fprintf(stderr,"Half edges full overlap:\n");
+						fprintf(stderr,"i = %d, j = %d\n",ep0->ix,ep1->ix);
+						fprintf(stderr,"%s %d at %c = %f from %c = %f to %f, half %s\n",
+						h == 0 ? "Vert" : "Horiz", ep0->ix,
+						h == 0 ? 'X' : 'Y', ep0->mj,
+						h == 0 ? 'Y' : 'X', ep0->mi0, ep0->mi1,
+						ep0->negh ? "Neg" : "Pos");
+						fprintf(stderr,"%s %d at %c = %f from %c = %f to %f, half %s\n",
+						h == 0 ? "Vert" : "Horiz", ep1->ix,
+						h == 0 ? 'X' : 'Y', ep1->mj,
+						h == 0 ? 'Y' : 'X', ep1->mi0, ep1->mi1,
+						ep1->negh ? "Neg" : "Pos");
+#endif
+						/* Split up the first edge so that we have four edges */
+						
+						if ((ep0b = (hedge *)calloc(sizeof(hedge), 1)) == NULL)
+							error("Malloc of half edge structure failed");
+						memcpy(ep0b, ep0, sizeof(hedge));
+
+						if ((ep0c = (hedge *)calloc(sizeof(hedge), 1)) == NULL)
+							error("Malloc of half edge structure failed");
+						memcpy(ep0c, ep0, sizeof(hedge));
+
+						ep0b->mi0 = ep1->mi0;
+						ep0b->mi1 = ep1->mi1;
+						ep0c->mi0 = ep1->mi1;
+						ep0->mi1 = ep1->mi0;
+
+						/* Insert them in order into linked list */
+						ep0c->next = ep1->next;
+						ep1->next = ep0c;
+						ep0b->next = ep1;
+						ep0->next = ep0b;
+						et.nhe[h] += 2;
+					}
+				}
+			}
+		}
+
+		/* Now we can assume that edges match completely or not at all */
+
+		/* Re-create sorted list of vertical half edges */
+		/* (Instead of this we could convert the half edge matching loops */
+		/*  below to use the linked list, like the above code.) */
+
+		if ((et.she[h] = (hedge **)malloc(sizeof(patch*) * et.nhe[h])) == NULL)
+			error("Malloc of array of vertical halfedge pointers failed");
+	
+		for (ep0 = et.he[h], i = 0; i < et.nhe[h]; i++, ep0 = ep0->next)
+			et.she[h][i] = ep0;
+	
+		/* Sort helf edges by their X location, then their Y0 location */
+#define HEAP_COMPARE(A,B) (fabs(A->mj - B->mj) < 1e-6 ? A->mi0 < B->mi0 : A->mj < B->mj)
+		HEAPSORT(hedge *, et.she[h], et.nhe[h]);
+#undef HEAP_COMPARE
 
 		et.nel[h] = 0;
 		et.nelp = &et.el[h];		/* Append next edge list here */
@@ -3980,7 +4143,7 @@ h == 0 ? 'Y' : 'X', et.she[h][i]->mi0, et.she[h][i]->mi1);
 			double *rgb = NULL;	/* Contrast RGB */
 			elist *el;		/* Current elist */
 
-			el = *et.nelp;
+			el = *et.nelp;	/* current end of lits */
 
 			/* Locate the end of the half edges at the same position */
 			for (ii = i; ii < et.nhe[h]; ii++) {
@@ -3988,11 +4151,13 @@ h == 0 ? 'Y' : 'X', et.she[h][i]->mi0, et.she[h][i]->mi1);
 					break;
 			}
 
-//printf("~1 doing group from %d to %d\n",i, ii);
+#ifdef DEBUG
+			fprintf(stderr,"Doing group from %d to %d\n",i, ii);
+#endif
+
 			/* Find half edge pairs */
 			/* Note that we assume that the half edges match perfectly, */
-			/* or not at all. This will be normaly be the case with targets */
-			/* generated by printtarg. */
+			/* or not at all. */
 			for (j = i; j < ii; j = nj, j++) {
 				int e, k = j+1;
 				double vv;
@@ -4013,25 +4178,26 @@ h == 0 ? 'Y' : 'X', et.she[h][i]->mi0, et.she[h][i]->mi1);
 					/* Found an overlapping non-matching edge */
 					nj = k;
 
-#ifdef NEVER
-fprintf(stderr,"i = %d, j = %d\n",i,j);
-fprintf(stderr,"%s %d at %c = %f from %c = %f to %f, half %s\n",
-h == 0 ? "Vert" : "Horiz", i,
-h == 0 ? 'X' : 'Y', et.she[h][j]->mj,
-h == 0 ? 'Y' : 'X', et.she[h][j]->mi0, et.she[h][j]->mi1,
-et.she[h][j]->negh ? "Neg" : "Pos");
-fprintf(stderr,"%s %d at %c = %f from %c = %f to %f, half %s\n",
-h == 0 ? "Vert" : "Horiz", i,
-h == 0 ? 'X' : 'Y', et.she[h][k]->mj,
-h == 0 ? 'Y' : 'X', et.she[h][k]->mi0, et.she[h][k]->mi1,
-et.she[h][k]->negh ? "Neg" : "Pos");
-#endif /* NEVER */
+#ifdef DEBUG
+					fprintf(stderr,"Half edges overlap but don't match:\n");
+					fprintf(stderr,"i = %d, j = %d\n",i,j);
+					fprintf(stderr,"%s %d at %c = %f from %c = %f to %f, half %s\n",
+					h == 0 ? "Vert" : "Horiz", i,
+					h == 0 ? 'X' : 'Y', et.she[h][j]->mj,
+					h == 0 ? 'Y' : 'X', et.she[h][j]->mi0, et.she[h][j]->mi1,
+					et.she[h][j]->negh ? "Neg" : "Pos");
+					fprintf(stderr,"%s %d at %c = %f from %c = %f to %f, half %s\n",
+					h == 0 ? "Vert" : "Horiz", j,
+					h == 0 ? 'X' : 'Y', et.she[h][k]->mj,
+					h == 0 ? 'Y' : 'X', et.she[h][k]->mi0, et.she[h][k]->mi1,
+					et.she[h][k]->negh ? "Neg" : "Pos");
+#endif /* DEBUG */
 						error("Internal - half edges don't match");
 
 				} else {
-					/* Must be a non-matching edge */
+					/* Must be a non-matching edge against the media */
 					nj = j;
-					rgb = et.mrgb;	/* Edge must be against media */
+					rgb = et.mrgb;
 				}
 
 				/* Compute vector delta in rgb */
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;
+}
+
 
 
 
diff --git a/target/targen.h b/target/targen.h
index 2f4665a..c99b551 100644
--- a/target/targen.h
+++ b/target/targen.h
@@ -23,6 +23,7 @@
 struct _fxpos {
 	double p[MXTD];		/* Device coordinate position */
 	double v[MXTD];		/* Room for perceptual value */
+	int eloc;			/* if >= 0, cgats index, to even location */
 }; typedef struct _fxpos fxpos;