Imported Upstream version 1.6.3upstream/1.6.3

1 files changed, 449 insertions, 176 deletions

diff --git a/spectro/dispcal.c b/spectro/dispcal.c
index accb79a..b8f3ebe 100644
--- a/spectro/dispcal.c
+++ b/spectro/dispcal.c
@@ -20,6 +20,21 @@
 
 /* TTBD
 
+	Calibrating the black point of a true power response
+	device is very slow to converge - the jacobian is always
+	underestimating the actual delta RGB needed because the
+	slope is getting shallower and shallower. Need to
+	be able to figure when to increase rgain in those circumstances,
+	rather than reducing it ? 
+
+	Dealing with noisy/inconsistent readings could probably
+	be improved - the statistical information from a iteration
+	series is being ignored. ie. do a linear regression/fit
+	on all the values for a given target, and then
+	at the end, use a weighted blend of the best solution
+	and the fit. Weight by something like the number used
+	for the fit. vs. 1.
+
 	Try to improve calibration speed by using adaptive
 	measurement set, rather than fixed resolution doubling ?
 	(ie. just measure at troublesome points using a "divide in half"
@@ -58,7 +73,7 @@
 	  display gives about 18% output at 50% device input.]
 
 
-	The verify (-E) may not be being done correctly.
+	The verify (-z) may not be being done correctly.
 	Like update, shouldn't it read the .cal file to set what's
 	being calibrated aganist ? (This would fix missing ambient value too!)
 
@@ -67,6 +82,13 @@
 	in "Figure out the black point target" - Yes they are !!
 	Verify probably shouldn't work this way.
 
+	Add DICOM support:
+
+		* Add 20% grey background full screen option + 10% patch recommendation
+		* Add "include Glare" option for contact instruments to dispsup.
+		* Add absolute DICOM function target.
+		* Add DICOM mode black point hue handling (? what policy ?)
+		* Add DICOM stats report (JND dE + mean + SD) to verify ??
  */
 
 #ifdef __MINGW32__
@@ -119,17 +141,22 @@
 
 #define COMPORT 1			/* Default com port 1..4 */
 #define OPTIMIZE_MODEL		/* Adjust model for best fit */
-#define REFINE_GAIN 0.80	/* Refinement correction damping/gain */
-#define MAX_RPTS 12			/* Maximum tries at making a sample meet the current threshold */
+#define REFINE_GAIN 0.90	/* Refinement correction damping/gain */
 #define VER_RES 100			/* Verification resolution */
 #define NEUTRAL_BLEND_RATE 4.0		/* Default rate of transition for -k factor < 1.0 (power) */
 #define ADJ_JACOBIAN		/* Adjust the Jacobian predictor matrix each time */
-#define JAC_COR_FACT 0.4	/* Amount to correct Jacobian by (to filter noise) */
-#define REMEAS_JACOBIAN		/* Re-measure Jacobian */
+#define JAC_COMP_FACT 0.4	/* Amount to compound Jacobian correction */
+#define JAC_COR_FACT 0.4	/* Amount to damp Jacobian by (to filter noise) */
+#define REMEAS_JACOBIAN		/* Re-measure Jacobian if it is a poor predictor */
 #define MOD_DIST_POW 1.6	/* Power used to distribute test samples for model building */
 #define REFN_DIST_POW 1.6	/* Power used to distribute test samples for grey axis refinement */
 #define CHECK_DIST_POW 1.6	/* Power used to distribute test samples for grey axis checking */
 #define THRESH_SCALE_POW 0.5 /* Amount to loosen threshold for first itterations */
+#define ADJ_THRESH			/* Adjust threshold to be half a step */
+#define MIN_THRESH 0.05	 	/* Minimum stopping threshold to allow in ADJ_THRESH */
+#define POWERR_THR 0.05		/* Point near black to start weighting +ve error */
+#define POWERR_WEIGHT 999.0	/* Weight to give +ve delta E at black */
+#define POWERR_WEIGHT_POW 4.0	/* Curve to plend from equal weight to +ve extra weight */
 #define CAL_RES 256			/* Resolution of calibration table to produce. */
 #define CLIP				/* Clip RGB during refinement */
 #define RDAC_SMOOTH 0.3		/* RAMDAC curve fitting smoothness */
@@ -190,6 +217,7 @@ typedef struct {
 
 	double nwh[3];		/* Target white normalised XYZ value (Y = 1.0) */
 	double twh[3];		/* Target white absolute XYZ value */
+	double twYxy[3];	/* Target white Yxy (informational) */
 	icmXYZNumber twN;	/* Same as above as XYZNumber */
 
 	double tbk[3];		/* Target black point color */
@@ -847,14 +875,16 @@ typedef struct {
 	double tXYZ[3];		/* Target XYZ */
 	double XYZ[3];		/* Read XYZ */
 	double deXYZ[3];	/* Delta XYZ wanted to target */
-	double de;			/* Delta Lab to neutral target */
-	double dc;			/* Delta XYZ to neutral target */
-	double peqde;		/* Delta Lab to previous point value (last pass) */
-	double hde;			/* Hybrid de composed of de and peqde */
+	double _de;			/* Non-weighted Delta Lab */
+	double de;			/* Weightd Delta Lab to neutral target */
+	double dc;			/* Weightd Delta XYZ to neutral target */
+	double peqde;		/* Weightd Delta Lab to last pass equivalent point value */
+	double hde;			/* Weightd Hybrid de composed of de and peqde */
 
+	double prgb[3];		/* Previous measured RGB */
 	double pXYZ[3];		/* Previous measured XYZ */
 	double pdXYZ[3];	/* Delta XYZ intended from previous measure */
-	double pdrgb[3];	/* Delta rgb made to previous */
+	double pdrgb[3];	/* Delta rgb made to previous to acorrect XYZ */
 
 	double dXYZ[3];		/* Actual delta XYZ resulting from previous delta rgb */
 
@@ -863,11 +893,12 @@ typedef struct {
 	double fb_ij[3][3];	/* Copy of initial inverse Jacobian, used as a fallback */
 } csp;
 
+
 /* All the sample points */
 typedef struct {
 	int no;				/* Number of samples */
 	int _no;			/* Allocation */
-	csp *s;			/* List of samples */
+	csp *s;				/* List of samples */
 } csamp;
 
 static void free_alloc_csamp(csamp *p) {
@@ -916,10 +947,13 @@ static void init_csamp_v(csamp *p, calx *x, int psrand) {
 }
 
 /* Initialise txyz values from v values */
-static void init_csamp_txyz(csamp *p, calx *x, int fixdev) {
+static void init_csamp_txyz(csamp *p, calx *x, int fixdev, int verb) {
 	int i, j;
 	double tbL[3];		/* tbk as Lab */
 
+	if (verb >= 3)
+		printf("init_csamp_txyz:\n");
+
 	/* Convert target black from XYZ to Lab here, */
 	/* in case twN has changed at some point. */
 	icmXYZ2Lab(&x->twN, tbL, x->tbk);
@@ -946,22 +980,22 @@ static void init_csamp_txyz(csamp *p, calx *x, int fixdev) {
 
 		/* Compute blended neutral target a* b* */
 		bl = pow((1.0 - vv), x->nbrate);		/* Crossover near the black */
-		Lab[1] = bl * tbL[1];
-		Lab[2] = bl * tbL[2];
+		Lab[1] = (1.0 - bl) * 0.0 + bl * tbL[1];
+		Lab[2] = (1.0 - bl) * 0.0 + bl * tbL[2];
 
 		icmAry2Ary(XYZ, p->s[i].tXYZ);				/* Save the existing values */
 		icmLab2XYZ(&x->twN, p->s[i].tXYZ, Lab);		/* New XYZ Value to aim for */
 
-#ifdef DEBUG
-		printf("%d: target XYZ %.2f %.2f %.2f, Lab %.2f %.2f %.2f\n",i, p->s[i].tXYZ[0],p->s[i].tXYZ[1],p->s[i].tXYZ[2], Lab[0],Lab[1],Lab[2]);
-#endif
+		if (verb >= 3) {
+			printf("%d: target XYZ %.4f %.4f %.4f, Lab %.3f %.3f %.3f\n",i, p->s[i].tXYZ[0],p->s[i].tXYZ[1],p->s[i].tXYZ[2], Lab[0],Lab[1],Lab[2]);
+		}
 	}
 }
 
 
 /* Allocate the sample points and initialise them with the */
 /* target device and XYZ values, and first cut device values. */
-static void init_csamp(csamp *p, calx *x, int doupdate, int verify, int psrand, int no) {
+static void init_csamp(csamp *p, calx *x, int doupdate, int verify, int psrand, int no, int verb) {
 	int i, j;
 	
 	p->_no = p->no = no;
@@ -971,7 +1005,7 @@ static void init_csamp(csamp *p, calx *x, int doupdate, int verify, int psrand,
 
 	/* Compute v and txyz */
 	init_csamp_v(p, x, psrand);
-	init_csamp_txyz(p, x, 0);
+	init_csamp_txyz(p, x, 0, verb);
 
 	/* Generate the sample points */
 	for (i = 0; i < no; i++) {
@@ -1064,7 +1098,7 @@ static void csamp_interp(csamp *p, double xyz[3], double v) {
 /* Re-initialise a CSP with a new number of points. */
 /* Interpolate the device values and jacobian. */
 /* Set the current rgb from the current RAMDAC curves if not verifying */
-static void reinit_csamp(csamp *p, calx *x, int verify, int psrand, int no) {
+static void reinit_csamp(csamp *p, calx *x, int verify, int psrand, int no, int verb) {
 	csp *os;			/* Old list of samples */
 	int ono;			/* Old number of samples */
 	int i, j, k, m;
@@ -1075,7 +1109,7 @@ static void reinit_csamp(csamp *p, calx *x, int verify, int psrand, int no) {
 	os = p->s;			/* Save the existing per point information */
 	ono = p->no;
 
-	init_csamp(p, x, 0, 2, psrand, no);
+	init_csamp(p, x, 0, 2, psrand, no, verb);
 
 	p->_no = p->no = no;
 
@@ -1151,7 +1185,7 @@ static void reinit_csamp(csamp *p, calx *x, int verify, int psrand, int no) {
 		/* Compute expected delta XYZ using new Jacobian */
 		icmMulBy3x3(p->s[i].pdXYZ, p->s[i].j, p->s[i].pdrgb);
 
-		p->s[i].de = b * os[j+1].de + (1.0 - b) * os[j].de;
+		p->s[i]._de = p->s[i].de = b * os[j+1].de + (1.0 - b) * os[j].de;
 		p->s[i].dc = b * os[j+1].dc + (1.0 - b) * os[j].dc;
 	}
 
@@ -1223,6 +1257,20 @@ static double comp_ct(
 	
 /* =================================================================== */
 
+/* Return the normal Delta E given two XYZ values, but */
+/* exagerate the L* error if act L* > targ L* by a factor of fact */
+extern ICCLIB_API double bwXYZLabDE(icmXYZNumber *w, double *targ, double *act, double fact) {
+	double targlab[3], actlab[3], rv;
+
+	icmXYZ2Lab(w, targlab, targ);
+	icmXYZ2Lab(w, actlab, act);
+	if (actlab[0] > targlab[0])
+		actlab[0] = targlab[0] + fact * (actlab[0] - targlab[0]);
+	rv = icmLabDE(targlab, actlab);
+	return rv;
+}
+/* =================================================================== */
+
 /* Default gamma */
 double g_def_gamma = 2.4;
 
@@ -1232,7 +1280,7 @@ double g_def_gamma = 2.4;
 
          ABCDEFGHIJKLMNOPQRSTUVWXYZ
   upper  .......... ....... . ...  
-  lower  .......   . ...... .... . 
+  lower  .......   . ...... .... ..
 
 */
 
@@ -1276,6 +1324,9 @@ void usage(char *diag, ...) {
 	}
 	free_disppaths(dp);
 	fprintf(stderr," -dweb[:port]         Display via a web server at port (default 8080)\n");
+#ifdef NT
+	fprintf(stderr," -dmadvr              Display via MadVR Video Renderer\n");
+#endif
 //	fprintf(stderr," -d fake              Use a fake display device for testing, fake%s if present\n",ICC_FILE_EXT);
 	fprintf(stderr," -c listno            Set communication port from the following list (default %d)\n",COMPORT);
 	if ((icmps = new_icompaths(g_log)) != NULL) {
@@ -1323,7 +1374,7 @@ void usage(char *diag, ...) {
 	fprintf(stderr," -A rate              Rate of blending from neutral to black point. Default %.1f\n",NEUTRAL_BLEND_RATE);
 	fprintf(stderr," -B blkbright         Set the target black brightness in cd/m^2\n");
 	fprintf(stderr," -e [n]               Run n verify passes on final curves\n");
-	fprintf(stderr," -E                   Run only verify pass on installed calibration curves\n");
+	fprintf(stderr," -z                   Run only verify pass on installed calibration curves\n");
 	fprintf(stderr," -P ho,vo,ss[,vs]     Position test window and scale it\n");
 	fprintf(stderr,"                      ho,vi: 0.0 = left/top, 0.5 = center, 1.0 = right/bottom etc.\n");
 	fprintf(stderr,"                      ss: 0.5 = half, 1.0 = normal, 2.0 = double etc.\n");
@@ -1331,6 +1382,7 @@ void usage(char *diag, ...) {
 #if defined(UNIX_X11)
 	fprintf(stderr," -n                   Don't set override redirect on test window\n");
 #endif
+	fprintf(stderr," -E                   Encode the test values for video range 16..235/255\n");
 	fprintf(stderr," -J                   Run instrument calibration first (used rarely)\n");
 	fprintf(stderr," -N                   Disable initial calibration of instrument if possible\n");
 	fprintf(stderr," -H                   Use high resolution spectrum mode (if available)\n");
@@ -1343,7 +1395,9 @@ void usage(char *diag, ...) {
 		fprintf(stderr,"                      1931_2 (def), 1964_10, S&B 1955_2, shaw, J&V 1978_2, 1964_10c\n");
 	}
 	fprintf(stderr," -I b|w               Drift compensation, Black: -Ib, White: -Iw, Both: -Ibw\n");
+	fprintf(stderr," -Y R:rate            Override measured refresh rate with rate Hz\n");
 	fprintf(stderr," -Y A                 Use non-adaptive integration time mode (if available).\n");
+	fprintf(stderr," -Y p                 Don't wait for the instrument to be placed on the display\n");
 	fprintf(stderr," -C \"command\"         Invoke shell \"command\" each time a color is set\n");
 	fprintf(stderr," -M \"command\"         Invoke shell \"command\" each time a color is measured\n");
 	fprintf(stderr," -W n|h|x             Override serial port flow control: n = none, h = HW, x = Xon/Xoff\n");
@@ -1360,6 +1414,7 @@ int main(int argc, char *argv[]) {
 	disppath *disp = NULL;				/* Display being used */
 	double hpatscale = 1.0, vpatscale = 1.0;	/* scale factor for test patch size */
 	double ho = 0.0, vo = 0.0;			/* Test window offsets, -1.0 to 1.0 */
+	int out_tvenc = 0;					/* 1 to use RGB Video Level encoding */
 	int blackbg = 0;            		/* NZ if whole screen should be filled with black */
 	int verb = 0;
 	int debug = 0;
@@ -1382,7 +1437,9 @@ int main(int argc, char *argv[]) {
 	int dtype = 0;						/* Display type selection charater */
 	int tele  = 0;						/* nz if telephoto mode */
 	int nocal = 0;						/* Disable auto calibration */
+	int noplace = 0;					/* Disable initial user placement check */
 	int highres = 0;					/* Use high res mode if available */
+	double refrate = 0.0;				/* 0.0 = default, > 0.0 = override refresh rate */ 
 	int nadaptive = 0;					/* Use non-adaptive mode if available */
 	int bdrift = 0;						/* Flag, nz for black drift compensation */
 	int wdrift = 0;						/* Flag, nz for white drift compensation */
@@ -1393,7 +1450,7 @@ int main(int argc, char *argv[]) {
 	double tbright = 0.0;				/* Target white brightness ( 0.0 == max)  */
 	double gamma = 0.0;					/* Advertised Gamma target */
 	double egamma = 0.0;				/* Effective Gamma target, NZ if set */
-	double ambient = 0.0;				/* NZ if viewing cond. adjustment to be used (cd/m^2) */
+	double ambient = 0.0;				/* NZ if viewing cond. adjustment to be used (Lux) */
 	double bkcorrect = -1.0;			/* Level of black point correction, < 0 = auto */ 
 	double bkbright = 0.0;				/* Target black brightness ( 0.0 == min)  */
 	int quality = -99;					/* Quality level, -2 = v, -1 = l, 0 = m, 1 = h, 2 = u */
@@ -1403,9 +1460,13 @@ int main(int argc, char *argv[]) {
 	int thrfail = 0;					/* Set to NZ if failed to meet threshold target */
 	double failerr = 0.0;				/* Delta E of worst failed target */
 	int mxits = 3;						/* maximum iterations (medium) */
+	int mxrpts = 12;					/* maximum repeats (medium) */
 	int verify = 0;						/* Do a verify after last refinement, 2 = do only verify. */
 	int nver = 0;						/* Number of verify passes after refinement */
 	int webdisp = 0;					/* NZ for web display, == port number */
+#ifdef NT
+	int madvrdisp = 0;					/* NZ for madvr display */
+#endif
 	char *ccallout = NULL;				/* Change color Shell callout */
 	char *mcallout = NULL;				/* Measure color Shell callout */
 	char outname[MAXNAMEL+1] = { 0 };	/* Output cgats file base name */
@@ -1421,9 +1482,12 @@ int main(int argc, char *argv[]) {
 	int it;								/* verify & refine iteration */
 	int rv;
 	int fitord = 30;					/* More seems to make curves smoother */
-	int native = 1;						/* 0 = use current or given calibration curve */
-										/* 1 = set native linear op and use ramdac high prec'n */
+	int native = 3;						/* X0 = use current per channel calibration curve */
+										/* X1 = set native linear output and use ramdac high prec */
+										/* 0X = use current color management cLut (MadVR) */
+										/* 1X = disable color management cLUT (MadVR) */
 	int noramdac = 0;					/* Will be set to nz if can't set ramdac */
+	int nocm = 0;						/* Will be set to nz if can't set color managament */
 	int errc;							/* Return value from new_disprd() */
 	calx x;								/* Context for calibration solution */
 
@@ -1513,6 +1577,12 @@ int main(int argc, char *argv[]) {
 							usage("Web port number must be in range 1..65535");
 					}
 					fa = nfa;
+#ifdef NT
+				} else if (strncmp(na,"madvr",5) == 0
+				 || strncmp(na,"MADVR",5) == 0) {
+					madvrdisp = 1;
+					fa = nfa;
+#endif
 				} else {
 #if defined(UNIX_X11)
 					int ix, iv;
@@ -1554,6 +1624,9 @@ int main(int argc, char *argv[]) {
 #endif
 				}
 
+			} else if (argv[fa][1] == 'E') {
+				out_tvenc = 1;
+
 			} else if (argv[fa][1] == 'J') {
 				docalib = 1;
 
@@ -1673,7 +1746,7 @@ int main(int argc, char *argv[]) {
 					fa = nfa;
 				}
 
-			} else if (argv[fa][1] == 'E') {
+			} else if (argv[fa][1] == 'z') {
 				verify = 2;
 				mfa = 0;
 
@@ -1834,7 +1907,6 @@ int main(int argc, char *argv[]) {
 				ambient = atof(na);
 				if (ambient < 0.0)
 					usage("-a parameter %f out of range",ambient);
-				ambient /= 3.141592654;	/* Convert from Lux to cd/m^2 */
 
 			/* Test patch offset and size */
 			} else if (argv[fa][1] == 'P') {
@@ -1864,11 +1936,20 @@ int main(int argc, char *argv[]) {
 				if (na == NULL)
 					usage("Flag '-Y' expects extra flag");
 			
-				if (na[0] == 'A') {
+				if (na[0] == 'R') {
+					if (na[1] != ':')
+						usage("-Y R:rate syntax incorrect");
+					refrate = atof(na+2);
+					if (refrate < 5.0 || refrate > 150.0)
+						usage("-Y R:rate %f Hz not in valid range",refrate);
+				} else if (na[0] == 'A') {
 					nadaptive = 1;
+				} else if (na[0] == 'p') {
+					noplace = 1;
 				} else {
 					usage("Flag '-Y %c' not recognised",na[0]);
 				}
+				fa = nfa;
 
 			} else 
 				usage("Flag '-%c' not recognised",argv[fa][1]);
@@ -1932,7 +2013,7 @@ int main(int argc, char *argv[]) {
 	}
 
 	if (fake)
-		comport = -99;
+		comport = FAKE_DEVICE_PORT;
 	if ((icmps = new_icompaths(g_log)) == NULL)
 		error("Finding instrument paths failed");
 	if ((ipath = icmps->get_path(icmps, comport)) == NULL)
@@ -1940,7 +2021,11 @@ int main(int argc, char *argv[]) {
 
 	if (docalib) {
 		if ((rv = disprd_calibration(ipath, fc, dtype, 0, tele, nadaptive, nocal, disp,
-		                             webdisp, blackbg, override,
+		                             webdisp,
+#ifdef NT
+			                         madvrdisp,
+#endif
+			                         out_tvenc, blackbg, override,
 			                         100.0 * hpatscale, 100.0 * vpatscale,
 			                         ho, vo, g_log)) != 0) {
 			error("docalibration failed with return value %d\n",rv);
@@ -1973,12 +2058,16 @@ int main(int argc, char *argv[]) {
 
 	/* Normally calibrate against native response */
 	if (verify == 2 || doreport == 2)
-		native = 0;	/* But measure calibrated response of verify or report calibrated */ 
+		native = 0;	/* But measure current calibrated & CM response for verify or report calibrated */ 
 
 	/* Get ready to do some readings */
-	if ((dr = new_disprd(&errc, ipath, fc, dtype, 0, tele, nadaptive, nocal,
-	                     highres, native, &noramdac, NULL, 0, 0, disp, blackbg, override,
-	                     webdisp, ccallout, mcallout,
+	if ((dr = new_disprd(&errc, ipath, fc, dtype, 0, tele, nadaptive, nocal, noplace,
+	                     highres, refrate, native, &noramdac, &nocm, NULL, 0,
+		                 disp, out_tvenc, blackbg, override, webdisp,
+#ifdef NT
+		                 madvrdisp,
+#endif
+		                 ccallout, mcallout,
 	                     100.0 * hpatscale, 100.0 * vpatscale, ho, vo,
 	                     cmx != NULL ? cmx->matrix : NULL,
 	                     ccs != NULL ? ccs->samples : NULL, ccs != NULL ? ccs->no_samp : 0,
@@ -1986,11 +2075,11 @@ int main(int argc, char *argv[]) {
 	                     "fake" ICC_FILE_EXT, g_log)) == NULL)
 		error("new_disprd() failed with '%s'\n",disprd_err(errc));
 
-	if (native != 0 && noramdac) {
-		warning("No access to VideoLUTs, so calibrating display as-is rather than native");
+	if ((native & 1) && noramdac) {
+		warning("Unable to access to VideoLUTs so can't be sure colors are native");
 		if (doprofile)
 			warning("Profile will reflect the as-is display response and not contain a 'vcgt' tag");
-		native = 0;
+		native &= ~1;
 
 		if (doupdate && doprofile)
 			error("Can't update a profile that doesn't use the 'vcgt' tag for calibration");
@@ -2053,7 +2142,7 @@ int main(int argc, char *argv[]) {
 		if (doreport == 1) {
 #define MAX_RES_SAMPS 24
 			col ttt[MAX_RES_SAMPS];
-			int res_samps = 6;
+			int res_samps = 9;
 			double a0, a1, a2, dd;
 			int n;
 			int issig = 0;
@@ -2073,7 +2162,7 @@ int main(int argc, char *argv[]) {
 					gcc_bug_fix(sigbits);
 #endif
 					/* Notional test value */
-					v = (5 << (sigbits-3))/((1 << sigbits) - 1.0);
+					v = (7 << (sigbits-3))/((1 << sigbits) - 1.0);
 					/* And -1, 0 , +1 bit test values */
 					if ((n % 3) == 2)
 						v += 1.0/((1 << sigbits) - 1.0);
@@ -2102,6 +2191,7 @@ int main(int argc, char *argv[]) {
 				a0 /= (res_samps / 3.0);
 				a1 /= (res_samps / 3.0);
 				a2 /= (res_samps / 3.0);
+				DBG((dbgo,"Bits %d: -1: %f 0: %f +1 %f\n",sigbits, a0, a1, a2));
 				/* Judge significance of any differences */
 				dd = 0.0;
 				for (n = 0; n < res_samps; n++) {
@@ -2120,13 +2210,13 @@ int main(int argc, char *argv[]) {
 					issig = 1;		/* Noticable difference */
 				else
 					issig = 0;		/* No noticable difference */
-				DBG((dbgo,"Bits %d: Between = %f, %f within = %f, sig = %s\n",sigbits, fabs(a1 - a0), fabs(a2 - a1),dd, issig ? "yes" : "no"));
+				DBG((dbgo,"Bits %d: Between = %f, %f within = %f, sig = %s\n",sigbits, fabs(a1 - a0), fabs(a2 - a1), dd, issig ? "yes" : "no"));
 
 				switch(sigbits) {
 					case 8:				/* Do another trial */
 						if (issig) {
 							sigbits = 10;
-							res_samps = 6;
+							res_samps = 9;
 						} else {
 							sigbits = 6;
 						}
@@ -2134,6 +2224,7 @@ int main(int argc, char *argv[]) {
 					case 6:				/* Do another trial or give up */
 						if (issig) {
 							sigbits = 7;
+							res_samps = 6;
 						} else {
 							sigbits = 0;
 							issig = 2;	/* Give up */
@@ -2147,7 +2238,7 @@ int main(int argc, char *argv[]) {
 					case 10:			/* Do another trial */
 						if (issig) {
 							sigbits = 12;
-							res_samps = 9;
+							res_samps = 12;
 						} else {
 							sigbits = 9;
 						}
@@ -2181,7 +2272,8 @@ int main(int argc, char *argv[]) {
 			printf("Current calibration response:\n");
 		else
 			printf("Uncalibrated response:\n");
-		printf("Black level = %.2f cd/m^2\n",tcols[0].XYZ[1]);
+		printf("Black level = %.4f cd/m^2\n",tcols[0].XYZ[1]);
+		printf("50%%   level = %.2f cd/m^2\n",tcols[1].XYZ[1]);
 		printf("White level = %.2f cd/m^2\n",tcols[2].XYZ[1]);
 		printf("Aprox. gamma = %.2f\n",cgamma);
 		printf("Contrast ratio = %.0f:1\n",tcols[2].XYZ[1]/tcols[0].XYZ[1]);
@@ -2239,6 +2331,13 @@ int main(int argc, char *argv[]) {
 			error("Can't update '%s' - there aren't two tables",outname);
 		}
 
+		out_tvenc = 0;
+		if ((fi = icg->find_kword(icg, 0, "TV_OUTPUT_ENCODING")) >= 0) {
+			if (strcmp(icg->t[0].kdata[fi], "YES") == 0
+			 || strcmp(icg->t[0].kdata[fi], "yes") == 0)
+				out_tvenc = 1;
+		}
+
 //printf("~1 reading previous cal, got 2 tables\n");
 
 		/* Read in the setup, user and model values */
@@ -2489,6 +2588,7 @@ int main(int argc, char *argv[]) {
 			isteps = 3;
 			rsteps = 9;
 			mxits = 1;
+			mxrpts = 8;
 			errthr = 2.0;
 			break;
 		case -2:				/* Very low */
@@ -2499,6 +2599,7 @@ int main(int argc, char *argv[]) {
 				mxits = 1;
 			else
 				mxits = 1;
+			mxrpts = 10;
 			break;
 		case -1:				/* Low */
 			if (verify != 2 && doprofile && !doupdate)
@@ -2511,6 +2612,7 @@ int main(int argc, char *argv[]) {
 				mxits = 1;
 			else
 				mxits = 2;
+			mxrpts = 10;
 			break;
 		default:
 		case 0:					/* Medum */
@@ -2525,6 +2627,7 @@ int main(int argc, char *argv[]) {
 				mxits = 1;
 			else
 				mxits = 3;
+			mxrpts = 12;
 			break;
 		case 1:					/* High */
 			if (verify != 2 && doprofile && !doupdate)
@@ -2537,6 +2640,7 @@ int main(int argc, char *argv[]) {
 				mxits = 1;
 			else
 				mxits = 4;
+			mxrpts = 16;
 			break;
 		case 2:					/* Ultra */
 			if (verify != 2 && doprofile && !doupdate)
@@ -2549,6 +2653,7 @@ int main(int argc, char *argv[]) {
 				mxits = 1;
 			else
 				mxits = 5;
+			mxrpts = 24;
 			break;
 	}
 
@@ -2560,6 +2665,9 @@ int main(int argc, char *argv[]) {
 
 	/* Say something about what we're doing */
 	if (verb) {
+		if (out_tvenc)
+			printf("Using TV encoding range of (16-235)/255\n");
+
 		if (dtype > 0)
 			printf("Display type is '%c'\n",dtype);
 
@@ -2665,11 +2773,11 @@ int main(int argc, char *argv[]) {
 				} 
 
 				if (verb) {
-					printf("Black = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0],
+					printf("Black = XYZ %6.4f %6.4f %6.4f\n",tcols[0].XYZ[0],
 					                                         tcols[0].XYZ[1], tcols[0].XYZ[2]);
-					printf("Grey  = XYZ %6.2f %6.2f %6.2f\n",tcols[1].XYZ[0],
+					printf("Grey  = XYZ %6.3f %6.3f %6.3f\n",tcols[1].XYZ[0],
 					                                         tcols[1].XYZ[1], tcols[1].XYZ[2]);
-					printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[2].XYZ[0],
+					printf("White = XYZ %6.3f %6.3f %6.3f\n",tcols[2].XYZ[0],
 					                                         tcols[2].XYZ[1], tcols[2].XYZ[2]);
 				}
 
@@ -2741,11 +2849,11 @@ int main(int argc, char *argv[]) {
 						error("display read failed with '%s'\n",disprd_err(rv));
 					}
 					if (verb) {
-						printf("Red   = XYZ %6.2f %6.2f %6.2f\n",ccols[0].XYZ[0],
+						printf("Red   = XYZ %6.3f %6.3f %6.3f\n",ccols[0].XYZ[0],
 						                       ccols[0].XYZ[1], ccols[0].XYZ[2]);
-						printf("Green = XYZ %6.2f %6.2f %6.2f\n",ccols[1].XYZ[0],
+						printf("Green = XYZ %6.3f %6.3f %6.3f\n",ccols[1].XYZ[0],
 						                       ccols[1].XYZ[1], ccols[1].XYZ[2]);
-						printf("Blue  = XYZ %6.2f %6.2f %6.2f\n",ccols[2].XYZ[0],
+						printf("Blue  = XYZ %6.3f %6.3f %6.3f\n",ccols[2].XYZ[0],
 						                       ccols[2].XYZ[1], ccols[2].XYZ[2]);
 					}
 					for (i = 0; i < 3; i++)
@@ -2758,7 +2866,7 @@ int main(int argc, char *argv[]) {
 					error("display read failed with '%s'\n",disprd_err(rv));
 				}
 				if (verb) {
-					printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0],
+					printf("White = XYZ %6.3f %6.3f %6.3f\n",tcols[0].XYZ[0],
 					                       tcols[0].XYZ[1], tcols[0].XYZ[2]);
 				}
 				
@@ -2939,7 +3047,7 @@ int main(int argc, char *argv[]) {
 					error("display read failed with '%s'\n",disprd_err(rv));
 				}
 				if (verb) {
-					printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0],
+					printf("White = XYZ %6.3f %6.3f %6.3f\n",tcols[0].XYZ[0],
 					                       tcols[0].XYZ[1], tcols[0].XYZ[2]);
 				}
 				
@@ -3012,11 +3120,11 @@ int main(int argc, char *argv[]) {
 						error("display read failed with '%s'\n",disprd_err(rv));
 					}
 					if (verb) {
-						printf("Red   = XYZ %6.2f %6.2f %6.2f\n",ccols[0].XYZ[0],
+						printf("Red   = XYZ %6.3f %6.3f %6.3f\n",ccols[0].XYZ[0],
 						                       ccols[0].XYZ[1], ccols[0].XYZ[2]);
-						printf("Green = XYZ %6.2f %6.2f %6.2f\n",ccols[1].XYZ[0],
+						printf("Green = XYZ %6.3f %6.3f %6.3f\n",ccols[1].XYZ[0],
 						                       ccols[1].XYZ[1], ccols[1].XYZ[2]);
-						printf("Blue  = XYZ %6.2f %6.2f %6.2f\n",ccols[2].XYZ[0],
+						printf("Blue  = XYZ %6.3f %6.3f %6.3f\n",ccols[2].XYZ[0],
 						                       ccols[2].XYZ[1], ccols[2].XYZ[2]);
 					}
 					for (i = 0; i < 3; i++)
@@ -3029,11 +3137,11 @@ int main(int argc, char *argv[]) {
 					error("display read failed with '%s'\n",disprd_err(rv));
 				}
 				if (verb) {
-					printf("Black = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0],
+					printf("Black = XYZ %6.4f %6.4f %6.4f\n",tcols[0].XYZ[0],
 					                                         tcols[0].XYZ[1], tcols[0].XYZ[2]);
-					printf("Grey  = XYZ %6.2f %6.2f %6.2f\n",tcols[1].XYZ[0],
+					printf("Grey  = XYZ %6.3f %6.3f %6.3f\n",tcols[1].XYZ[0],
 					                                         tcols[1].XYZ[1], tcols[1].XYZ[2]);
-					printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[2].XYZ[0],
+					printf("White = XYZ %6.3f %6.3f %6.3f\n",tcols[2].XYZ[0],
 					                                         tcols[2].XYZ[1], tcols[2].XYZ[2]);
 				}
 				
@@ -3064,7 +3172,7 @@ int main(int argc, char *argv[]) {
 				}
 
 				printf("\nAdjust R,G & B offsets to get target x,y. Press space when done.\n");
-				printf("   Target Br %.2f, x %.4f , y %.4f \n", tar1, tYxy[1],tYxy[2]);
+				printf("   Target Br %.4f, x %.4f , y %.4f \n", tar1, tYxy[1],tYxy[2]);
 				for (ff = 0;; ff ^= 1) {
 					double dir;			/* Direction to adjust brightness */
 					double sv[3], val1;				/* Scaled values */
@@ -3133,7 +3241,7 @@ int main(int argc, char *argv[]) {
 						rgbdir[0] = rgbdir[1] = rgbdir[2] = 0.0;
 					rgbxdir[bx] = rgbdir[bx];
 
-			 		printf("%c%c Current Br %.2f, x %.4f%c, y %.4f%c  DE %4.1f  R%c%c G%c%c B%c%c ",
+			 		printf("%c%c Current Br %.4f, x %.4f%c, y %.4f%c  DE %4.1f  R%c%c G%c%c B%c%c ",
 					       cr_char,
 					       ff == 0 ? '/' : '\\',
 					       val1,
@@ -3177,11 +3285,11 @@ int main(int argc, char *argv[]) {
 					error("display read failed with '%s'\n",disprd_err(rv));
 				}
 				if (verb) {
-					printf("Black = XYZ %6.2f %6.2f %6.2f\n",tcols[0].XYZ[0],
+					printf("Black = XYZ %6.4f %6.4f %6.4f\n",tcols[0].XYZ[0],
 					                                         tcols[0].XYZ[1], tcols[0].XYZ[2]);
-					printf("Grey  = XYZ %6.2f %6.2f %6.2f\n",tcols[1].XYZ[0],
+					printf("Grey  = XYZ %6.3f %6.3f %6.3f\n",tcols[1].XYZ[0],
 					                                         tcols[1].XYZ[1], tcols[1].XYZ[2]);
-					printf("White = XYZ %6.2f %6.2f %6.2f\n",tcols[2].XYZ[0],
+					printf("White = XYZ %6.3f %6.3f %6.3f\n",tcols[2].XYZ[0],
 					                                         tcols[2].XYZ[1], tcols[2].XYZ[2]);
 				}
 				
@@ -3198,7 +3306,7 @@ int main(int argc, char *argv[]) {
 					error("display read failed with '%s'\n",disprd_err(rv));
 				}
 				if (verb) {
-					printf("1%%    = XYZ %6.2f %6.2f %6.2f\n",tcols[3].XYZ[0],
+					printf("1%%    = XYZ %6.3f %6.3f %6.3f\n",tcols[3].XYZ[0],
 					                                         tcols[3].XYZ[1], tcols[3].XYZ[2]);
 				}
 				
@@ -3277,17 +3385,17 @@ int main(int argc, char *argv[]) {
 				printf("\n");
 
 				if (tbright > 0.0)			/* Given brightness */
-					printf("  Target Brightness = %.2f, Current = %5.2f, error = % .1f%%\n",
+					printf("  Target Brightness = %.3f, Current = %.3f, error = % .1f%%\n",
 				       tarw, tcols[2].XYZ[1], 
 				       100.0 * (tcols[2].XYZ[1] - tarw)/tarw);
 				else
 					printf("  Current Brightness = %.2f\n", tcols[2].XYZ[1]);
 				
-				printf("  Target 50%% Level  = %.2f, Current = %5.2f, error = % .1f%%\n",
+				printf("  Target 50%% Level  = %.3f, Current = %.3f, error = % .1f%%\n",
 				       tarh, tcols[1].XYZ[1], 
 				       100.0 * (tcols[1].XYZ[1] - tarh)/tarw);
 				
-				printf("  Target Near Black = %5.2f, Current = %5.2f, error = % .1f%%\n",
+				printf("  Target Near Black = %.4f, Current = %.4f, error = % .1f%%\n",
 				       tar1, val1, 
 				       100.0 * (val1 - tar1)/tarw);
 
@@ -3312,7 +3420,7 @@ int main(int argc, char *argv[]) {
 						error("ambient measure failed with '%s'\n",disprd_err(rv));
 					}
 				} else {
-					printf("Measured ambient level = %.1f Lux\n",ambient * 3.141592654);
+					printf("Measured ambient level = %.1f Lux\n",ambient);
 				}
 
 			} else if (c == '7') {
@@ -3349,14 +3457,26 @@ int main(int argc, char *argv[]) {
 		icmXYZNumber mwh;
 		ramdac *or = NULL;
 
-		col base[6] = {	/* Base set of test colors */
+		col base[9] = {	/* Base set of test colors */
 			{ 0.0, 0.0, 0.0 },		/* 0 - Black */
 			{ 1.0, 0.0, 0.0 },		/* 1 - Red */
-			{ 0.0, 1.0, 0.0 },		/* 2 - Green */
-			{ 0.0, 0.0, 1.0 },		/* 3 - Blue */
-			{ 1.0, 1.0, 1.0 },		/* 4 - White */
-			{ 0.0, 0.0, 0.0 }		/* 5 - Black */
+			{ 1.0, 1.0, 1.0 },		/* 2 - White */
+			{ 0.0, 0.0, 0.0 },		/* 3 - Black */
+			{ 0.0, 1.0, 0.0 },		/* 4 - Green */
+			{ 0.0, 0.0, 0.0 },		/* 5 - Black */
+			{ 0.0, 0.0, 1.0 },		/* 6 - Blue */
+			{ 1.0, 1.0, 1.0 },		/* 7 - White */
+			{ 0.0, 0.0, 0.0 }		/* 8 - Black */
 		};
+		int ix_k1 = 0;
+		int ix_r  = 1;
+		int ix_w1 = 2;
+		int ix_k2 = 3;
+		int ix_g  = 4;
+		int ix_k3 = 5;
+		int ix_b  = 6;
+		int ix_w2 = 7;
+		int ix_k4 = 8;
 
 		if (verb) {
 			if (verify == 2)
@@ -3379,7 +3499,7 @@ int main(int argc, char *argv[]) {
 		}
 
 		/* Read the patches without clamping */
-		if ((rv = dr->read(dr, base, 6, 1, 6, 1, 0, instNoClamp)) != 0) {
+		if ((rv = dr->read(dr, base, 9, 1, 9, 1, 0, instNoClamp)) != 0) {
 			dr->del(dr);
 			error("display read failed with '%s'\n",disprd_err(rv));
 		} 
@@ -3391,35 +3511,49 @@ int main(int argc, char *argv[]) {
 			or->del(or);
 		}
 
-		if (base[0].XYZ_v == 0) {
+		if (base[ix_k1].XYZ_v == 0) {
 			dr->del(dr);
 			error("Failed to get an XYZ value from the instrument!\n");
 		}
 
-		/* Average black relative from 2 readings */
-		x.bk[0] = 0.5 * (base[0].XYZ[0] + base[5].XYZ[0]);
-		x.bk[1] = 0.5 * (base[0].XYZ[1] + base[5].XYZ[1]);
-		x.bk[2] = 0.5 * (base[0].XYZ[2] + base[5].XYZ[2]);
+		if (verb >= 3) {
+			for (i = 0; i < 9; i++)
+				printf("Meas %d XYZ = %f %f %f\n",i,base[i].XYZ[0], base[i].XYZ[1], base[i].XYZ[2]);
+		}
+
+		/* Average black relative from 4 readings */
+		x.bk[0] = 0.25 * (base[ix_k1].XYZ[0] + base[ix_k2].XYZ[0]
+		                + base[ix_k3].XYZ[0] + base[ix_k4].XYZ[0]);
+		x.bk[1] = 0.25 * (base[ix_k1].XYZ[1] + base[ix_k2].XYZ[1]
+		                + base[ix_k3].XYZ[1] + base[ix_k4].XYZ[1]);
+		x.bk[2] = 0.25 * (base[ix_k1].XYZ[2] + base[ix_k2].XYZ[2]
+		                + base[ix_k3].XYZ[2] + base[ix_k4].XYZ[2]);
 		icmClamp3(x.bk, x.bk);	/* And clamp them */
-		for (i = 1; i < 5; i++)
+
+		/* Average white reading from 2 readings */
+		base[ix_w1].XYZ[0] = 0.5 * (base[ix_w1].XYZ[0] + base[ix_w2].XYZ[0]);
+		base[ix_w1].XYZ[1] = 0.5 * (base[ix_w1].XYZ[1] + base[ix_w2].XYZ[1]);
+		base[ix_w1].XYZ[2] = 0.5 * (base[ix_w1].XYZ[2] + base[ix_w2].XYZ[2]);
+	
+		for (i = 0; i < 9; i++)
 			icmClamp3(base[i].XYZ, base[i].XYZ);
 
 		/* Copy other readings into place */
-		dispLum = base[4].XYZ[1];				/* White Y */
-		icmAry2Ary(x.wh, base[4].XYZ);
+		dispLum = base[ix_w1].XYZ[1];				/* White Y */
+		icmAry2Ary(x.wh, base[ix_w1].XYZ);
 		icmAry2XYZ(x.twN, x.wh);	/* Use this as Lab reference white until we establish target */
 
-		icmAry2XYZ(mrd, base[1].XYZ);
-		icmAry2XYZ(mgn, base[2].XYZ);
-		icmAry2XYZ(mbl, base[3].XYZ);
-		icmAry2XYZ(mwh, base[4].XYZ);
+		icmAry2XYZ(mrd, base[ix_r].XYZ);
+		icmAry2XYZ(mgn, base[ix_g].XYZ);
+		icmAry2XYZ(mbl, base[ix_b].XYZ);
+		icmAry2XYZ(mwh, base[ix_w1].XYZ);
 
 		if (verb) {
-			printf("Black = XYZ %6.2f %6.2f %6.2f\n",x.bk[0],x.bk[1],x.bk[2]);
-			printf("Red   = XYZ %6.2f %6.2f %6.2f\n",base[1].XYZ[0], base[1].XYZ[1], base[1].XYZ[2]);
-			printf("Green = XYZ %6.2f %6.2f %6.2f\n",base[2].XYZ[0], base[2].XYZ[1], base[2].XYZ[2]);
-			printf("Blue  = XYZ %6.2f %6.2f %6.2f\n",base[3].XYZ[0], base[3].XYZ[1], base[3].XYZ[2]);
-			printf("White = XYZ %6.2f %6.2f %6.2f\n",base[4].XYZ[0], base[4].XYZ[1], base[4].XYZ[2]);
+			printf("Black = XYZ %6.4f %6.4f %6.4f\n",x.bk[0],x.bk[1],x.bk[2]);
+			printf("Red   = XYZ %6.3f %6.3f %6.3f\n",base[ix_r].XYZ[0], base[ix_r].XYZ[1], base[ix_r].XYZ[2]);
+			printf("Green = XYZ %6.3f %6.3f %6.3f\n",base[ix_g].XYZ[0], base[ix_g].XYZ[1], base[ix_g].XYZ[2]);
+			printf("Blue  = XYZ %6.3f %6.3f %6.3f\n",base[ix_b].XYZ[0], base[ix_b].XYZ[1], base[ix_b].XYZ[2]);
+			printf("White = XYZ %6.3f %6.3f %6.3f\n",base[ix_w1].XYZ[0], base[ix_w1].XYZ[1], base[ix_w1].XYZ[2]);
 		}
 
 		/* Setup forward matrix */
@@ -3427,6 +3561,7 @@ int main(int argc, char *argv[]) {
 			dr->del(dr);
 			error("Aprox. fwd matrix unexpectedly singular\n");
 		}
+		icmTranspose3x3(x.fm, x.fm);	/* Convert [RGB][XYZ] to [XYZ][RGB] */
 
 #ifdef DEBUG
 		if (verb) {
@@ -3783,8 +3918,11 @@ int main(int argc, char *argv[]) {
 			}
 		}
 
+		icmXYZ2Yxy(x.twYxy, x.twh);		/* For information */
+
 		if (verb)
-			printf("Target white value is XYZ %f %f %f\n",x.twh[0],x.twh[1],x.twh[2]);
+			printf("Target white value is XYZ %f %f %f [xy %f %f]\n",x.twh[0],x.twh[1],x.twh[2],
+			                                                          x.twYxy[1], x.twYxy[2]);
 	}
 
 	/* Need this for Lab conversions */
@@ -3800,6 +3938,7 @@ int main(int argc, char *argv[]) {
 //printf("~1 black point Lab = %f %f %f\n", tbkLab[0], tbkLab[1], tbkLab[2]);
 
 		/* Now blend the a* b* with that of the target white point */
+		/* according to how much to try and correct the hue. */
 		tbL[0] = tbkLab[0];
 		tbL[1] = bkcorrect * 0.0 + (1.0 - bkcorrect) * tbkLab[1];
 		tbL[2] = bkcorrect * 0.0 + (1.0 - bkcorrect) * tbkLab[2];
@@ -3829,7 +3968,7 @@ int main(int argc, char *argv[]) {
 		icmLab2XYZ(&x.twN, x.tbk, tbL);
 		icmAry2XYZ(x.tbN, x.tbk);
 		if (verb)
-			printf("Adjusted target black XYZ %.2f %.2f %.2f, Lab %.2f %.2f %.2f\n",
+			printf("Adjusted target black XYZ %.4f %.4f %.4f, Lab %.3f %.3f %.3f\n",
 	        x.tbk[0], x.tbk[1], x.tbk[2], tbL[0], tbL[1], tbL[2]);
 	}
 
@@ -3857,11 +3996,17 @@ int main(int argc, char *argv[]) {
 
 		if (verb) {
 			double tbp[3], tbplab[3];
-	        tbp[0] = x.tbk[0] * tby/x.tbk[1];
+			if (fabs(x.tbk[1]) > 1e-9)
+		        tbp[0] = x.tbk[0] * tby/x.tbk[1];
+			else
+		        tbp[0] = x.tbk[0];
 			tbp[1] = tby;
-			tbp[2] = x.tbk[2] * tby/x.tbk[1];
+			if (fabs(x.tbk[1]) > 1e-9)
+				tbp[2] = x.tbk[2] * tby/x.tbk[1];
+			else
+				tbp[2] = x.tbk[2];
 			icmXYZ2Lab(&x.twN, tbplab, tbp);
-			printf("Target black after min adjust: XYZ %.3f %.3f %.3f, Lab %.3f %.3f %.3f\n",
+			printf("Target black after min adjust: XYZ %.4f %.4f %.4f, Lab %.3f %.3f %.3f\n",
 			        tbp[0], tbp[1], tbp[2], tbplab[0], tbplab[1], tbplab[2]);
 		}
 
@@ -3896,7 +4041,7 @@ int main(int argc, char *argv[]) {
 					0.2 * 80.0,				/* Adapting luminence, 20% of display 80 cd/m^2 */
 					0.2,					/* Background relative to reference white */
 					80.0,					/* Display is 80 cd/m^2 */
-			        0.01, x.nwh,			/* 1% flare same white point */
+			        0.0, 0.01, x.nwh,		/* 0% flare and 1% glare same white point */
 					0);
 				break;
 
@@ -3907,7 +4052,7 @@ int main(int argc, char *argv[]) {
 					0.2 * 1000.0/3.1415,	/* Adapting luminence, 20% of 1000 lux in cd/m^2 */
 					0.2,					/* Background relative to reference white */
 					1000.0/3.1415,			/* Luminance of white in the Image field (cd/m^2) */
-			        0.01, x.nwh,			/* 1% flare same white point */
+			        0.0, 0.01, x.nwh,		/* 0% flare and 1% glare same white point */
 					0);
 				break;
 
@@ -3917,10 +4062,10 @@ int main(int argc, char *argv[]) {
 		/* The display we're calibratings situation */
 		x.dvc->set_view(x.dvc, vc_none,
 			x.nwh,				/* Display normalised white point */
-			0.2 * ambient,		/* Adapting luminence, 20% of ambient in cd/m^2 */
+			0.2 * ambient/3.1415,	/* Adapting luminence, 20% of ambient in cd/m^2 */
 			0.2,				/* Background relative to reference white */
 			x.twh[1],			/* Target white level (cd/m^2) */
-	        0.01, x.nwh,		/* 1% flare same white point */
+	        0.0, 0.01, x.nwh,	/* 0% flare and 1% glare same white point */
 			0);
 
 		/* Compute the normalisation values */
@@ -3984,7 +4129,7 @@ int main(int argc, char *argv[]) {
 	}
 
 	/* Setup the initial calibration test point values */
-	init_csamp(&asgrey, &x, doupdate, verify, verify == 2 ? 1 : 0, rsteps);
+	init_csamp(&asgrey, &x, doupdate, verify, verify == 2 ? 1 : 0, rsteps, verb);
 	
 	/* Calculate the initial calibration curve values */
 	if (verify != 2 && !doupdate) {
@@ -4069,13 +4214,13 @@ int main(int argc, char *argv[]) {
 
 		/* Verify pass */
 		if (it >= mxits)
-			rsteps = VER_RES;		/* Fixed verification resolution */
+			rsteps = VER_RES;	/* Fixed verification resolution */
 		else
-			thrfail = 0; /* Not verify pass */
+			thrfail = 0;		/* Not verify pass */
 
 
 		/* re-init asgrey if the number of test points has changed */
-		reinit_csamp(&asgrey, &x, verify, (verify == 2 || it >= mxits) ? 1 : 0, rsteps);
+		reinit_csamp(&asgrey, &x, verify, (verify == 2 || it >= mxits) ? 1 : 0, rsteps, verb);
 
 		if (verb) {
 			if (it >= mxits)
@@ -4084,41 +4229,67 @@ int main(int argc, char *argv[]) {
 				printf("Doing iteration %d with %d sample points and repeat threshold of %f DE\n",
 				                                                             it+1,rsteps, errthr);
 		}
+
 		/* Read and adjust each step */
 		/* Do this white to black to track drift in native white point */
 		for (i = rsteps-1; i >= 0; i--) {
-			double rpt;
+			int rpt;
 			double peqXYZ[3];		/* Previous steps equivalent aim point */
 			double bestrgb[3];		/* In case we fail */
 			double bestxyz[3];
 			double prevde = 1e7;
+			double best_de = 1e7;
 			double bestde = 1e7;
 			double bestdc = 1e7;
 			double bestpeqde = 1e7;
 			double besthde = 1e7;
 			double rgain = REFINE_GAIN;	/* Scale down if lots of repeats */
 			int mjac = 0;				/* We measured the Jacobian */
+			double ierrth = errthr;		/* This points error threshold */
+			
+			icmCpy3(asgrey.s[i].prgb, asgrey.s[i].rgb);		/* Init previous */
 
-			/* Setup a second termination threshold chriteria based on */
+			/* Setup a second termination threshold criteria based on */
 			/* the delta E to the previous step point for the last pass. */
+			/* This is to try and steer towards neutral axis consistency ? */
 			if (i == (rsteps-1) || it < (mxits-1)) {
 				icmAry2Ary(peqXYZ, asgrey.s[i].tXYZ);	/* Its own aim point */
 			} else {
 				double Lab1[3], Lab2[3], Lab3[3];
 				icmXYZ2Lab(&x.twN, Lab1, asgrey.s[i+1].XYZ);
 				icmXYZ2Lab(&x.twN, Lab2, asgrey.s[i].tXYZ);
-				Lab3[0] = Lab2[0];
-				Lab3[1] = 0.5 * (Lab1[1] + Lab2[1]);	/* Compute aim point between actual target */
-				Lab3[2] = 0.5 * (Lab1[2] + Lab2[2]);	/* and previous point. */
-				icmLab2XYZ(&x.twN, asgrey.s[i].tXYZ, Lab3);
 				Lab1[0] = Lab2[0];		/* L of current target with ab of previous as 2nd threshold */
-				icmLab2XYZ(&x.twN, peqXYZ, Lab1);
+				icmLab2XYZ(&x.twN, peqXYZ, Lab1);		/* Previous equivalent */
 			}
 
+#ifdef ADJ_THRESH
+			/* Adjust the termination threshold to make sure it is less than */
+			/* half a step */
+			{
+				double de;
+				if (i < (rsteps-1)) {
+					de = 0.5 * icmXYZLabDE(&x.twN, asgrey.s[i].tXYZ, asgrey.s[i+1].tXYZ);
+					if (de < MIN_THRESH)		/* Don't be silly */
+						de = MIN_THRESH;
+					if (de < ierrth)
+						ierrth = de;
+				}
+				if (i > 0) {
+					de = 0.5 * icmXYZLabDE(&x.twN, asgrey.s[i].tXYZ, asgrey.s[i-1].tXYZ);
+					if (de < MIN_THRESH)		/* Don't be silly */
+						de = MIN_THRESH;
+					if (de < ierrth)
+						ierrth = de;
+				}
+			}
+#endif /* ADJ_THRESH */
+
 			/* Until we meet the necessary accuracy or give up */
-			for (rpt = 0;rpt < MAX_RPTS; rpt++) {
+			for (rpt = 0; rpt < mxrpts; rpt++) {
+				double hlew = 1.0;	/* high L* error weight */
 				int gworse = 0;		/* information flag */
-				double wde;			/* informational */
+				double w_de, wde;	/* informational */
+				double pjadj[3][3] = { 0.0 };	/* Previous jacobian adjustment */
 
 				set[0].r = asgrey.s[i].rgb[0];
 				set[0].g = asgrey.s[i].rgb[1];
@@ -4142,7 +4313,7 @@ int main(int argc, char *argv[]) {
 
 					icmAry2Ary(x.twh, asgrey.s[i].XYZ);		/* Set current white */
 					icmAry2XYZ(x.twN, x.twh);				/* Need this for Lab conversions */
-					init_csamp_txyz(&asgrey, &x, 1);		/* Recompute txyz's */
+					init_csamp_txyz(&asgrey, &x, 1, verb);	/* Recompute txyz's */
 					icmAry2Ary(peqXYZ, asgrey.s[i].tXYZ);	/* Fix peqXYZ */
 //printf("~1 Just reset target white to native white\n");
 					if (wdrift) {	/* Make sure white drift is reset on next read. */
@@ -4150,35 +4321,48 @@ int main(int argc, char *argv[]) {
 					}
 				}
 
-				/* Compute the current change wanted to hit target */
+				/* Compute the next change wanted to hit target */
 				icmSub3(asgrey.s[i].deXYZ, asgrey.s[i].tXYZ, asgrey.s[i].XYZ);
-				asgrey.s[i].de = icmXYZLabDE(&x.twN, asgrey.s[i].tXYZ, asgrey.s[i].XYZ);
-				asgrey.s[i].peqde = icmXYZLabDE(&x.twN, peqXYZ, asgrey.s[i].XYZ);
+
+				/* For the darkest 5% of targets, weight the L* delta E so that */
+				/* we err on the darker side */
+				if (asgrey.s[i].v < POWERR_THR)
+					hlew = 1.0 + POWERR_WEIGHT * pow((POWERR_THR - asgrey.s[i].v)/POWERR_THR,
+					                                                       POWERR_WEIGHT_POW);
+				else
+					hlew = 1.0;
+//printf("~1 i %d, v %f, hlew %f\n",i,asgrey.s[i].v,hlew);
+				asgrey.s[i]._de = icmXYZLabDE(&x.twN, asgrey.s[i].tXYZ, asgrey.s[i].XYZ);
+				asgrey.s[i].de = bwXYZLabDE(&x.twN, asgrey.s[i].tXYZ, asgrey.s[i].XYZ, hlew);
+				asgrey.s[i].peqde = bwXYZLabDE(&x.twN, peqXYZ, asgrey.s[i].XYZ, hlew);
 				asgrey.s[i].hde = 0.8 * asgrey.s[i].de + 0.2 * asgrey.s[i].peqde;
 				/* Eudclidian difference of XYZ, because this doesn't always track Lab */
 				asgrey.s[i].dc = icmLabDE(asgrey.s[i].tXYZ, asgrey.s[i].XYZ);
 
-				/* Compute change from last XYZ */
+				/* Compute actual change from last XYZ */
 				icmSub3(asgrey.s[i].dXYZ, asgrey.s[i].XYZ, asgrey.s[i].pXYZ);
 
-#ifdef DEBUG
-				printf("\n\nTest point %d, v = %f\n",rsteps - i,asgrey.s[i].v);
-				printf("Current rgb %f %f %f -> XYZ %f %f %f, de %f, dc %f\n", 
-				asgrey.s[i].rgb[0], asgrey.s[i].rgb[1], asgrey.s[i].rgb[2],
-				asgrey.s[i].XYZ[0], asgrey.s[i].XYZ[1], asgrey.s[i].XYZ[2],
-				asgrey.s[i].de, asgrey.s[i].dc);
-				printf("Target XYZ %f %f %f, delta needed %f %f %f\n", 
-				asgrey.s[i].tXYZ[0], asgrey.s[i].tXYZ[1], asgrey.s[i].tXYZ[2],
-				asgrey.s[i].deXYZ[0], asgrey.s[i].deXYZ[1], asgrey.s[i].deXYZ[2]);
-				if (rpt > 0) {
-					printf("Intended XYZ change %f %f %f, actual change %f %f %f\n", 
-					asgrey.s[i].pdXYZ[0], asgrey.s[i].pdXYZ[1], asgrey.s[i].pdXYZ[2],
-					asgrey.s[i].dXYZ[0], asgrey.s[i].dXYZ[1], asgrey.s[i].dXYZ[2]);
+				w_de = asgrey.s[i]._de;
+				wde = asgrey.s[i].de;
+
+				if (verb >= 3) {
+					printf("\n\nTest point %d, v = %f, rpt %d\n",rsteps - i,asgrey.s[i].v,rpt);
+					printf("Current rgb %f %f %f -> XYZ %f %f %f, de %f, dc %f\n", 
+					asgrey.s[i].rgb[0], asgrey.s[i].rgb[1], asgrey.s[i].rgb[2],
+					asgrey.s[i].XYZ[0], asgrey.s[i].XYZ[1], asgrey.s[i].XYZ[2],
+					asgrey.s[i]._de, asgrey.s[i].dc);
+					printf("Target XYZ %f %f %f, delta needed %f %f %f\n", 
+					asgrey.s[i].tXYZ[0], asgrey.s[i].tXYZ[1], asgrey.s[i].tXYZ[2],
+					asgrey.s[i].deXYZ[0], asgrey.s[i].deXYZ[1], asgrey.s[i].deXYZ[2]);
+					if (rpt > 0) {
+						printf("Last intended XYZ change %f %f %f, actual change %f %f %f\n", 
+						asgrey.s[i].pdXYZ[0], asgrey.s[i].pdXYZ[1], asgrey.s[i].pdXYZ[2],
+						asgrey.s[i].dXYZ[0], asgrey.s[i].dXYZ[1], asgrey.s[i].dXYZ[2]);
+					}
 				}
-#endif
 
 				if (it < mxits) {		/* Not verify, apply correction */
-					int impj = 0;		/* We improved the Jacobian */
+					int impj = 0;		/* We adjusted the Jacobian */
 					int dclip = 0;		/* We clipped the new RGB */
 #ifdef ADJ_JACOBIAN
 					int isclipped = 0;
@@ -4197,7 +4381,9 @@ int main(int argc, char *argv[]) {
 
 #ifdef REMEAS_JACOBIAN
 					/* If the de hasn't improved, try and measure the Jacobian */
-					if (it < (rsteps-1) && mjac == 0 && asgrey.s[i].de > (0.8 * prevde)) {
+//					if (it < (rsteps-1) && mjac == 0 && asgrey.s[i].de > (0.8 * prevde))
+					if (mjac == 0 && asgrey.s[i].de > (0.8 * prevde))
+					{
 						double dd;
 						if (asgrey.s[i].v < 0.5)
 							dd = 0.05;
@@ -4224,18 +4410,20 @@ int main(int argc, char *argv[]) {
 						} 
 						totmeas += 3;
 
+//printf("\n~1 remeasured jacobian\n");
 						/* Matrix organization is J[XYZ][RGB] for del RGB->del XYZ*/
 						for (j = 0; j < 3; j++) {
 							asgrey.s[i].j[0][j] = (set[j].XYZ[0] - asgrey.s[i].XYZ[0]) / dd;
 							asgrey.s[i].j[1][j] = (set[j].XYZ[1] - asgrey.s[i].XYZ[1]) / dd;
 							asgrey.s[i].j[2][j] = (set[j].XYZ[2] - asgrey.s[i].XYZ[2]) / dd;
 						}
+
+						/* Clear pjadj */
+						for (j = 0; j < 3; j++)
+							pjadj[3][0] = pjadj[3][1] = pjadj[3][2] = 0.0;
+
 						if (icmInverse3x3(asgrey.s[i].ij, asgrey.s[i].j)) {
 							/* Should repeat with bigger dd ? */
-//printf("~1 matrix =\n");
-//printf("~1    %f %f %f\n", asgrey.s[i].j[0][0], asgrey.s[i].j[0][1], asgrey.s[i].j[0][2]);
-//printf("~1    %f %f %f\n", asgrey.s[i].j[1][0], asgrey.s[i].j[1][1], asgrey.s[i].j[1][2]);
-//printf("~1    %f %f %f\n", asgrey.s[i].j[2][0], asgrey.s[i].j[2][1], asgrey.s[i].j[2][2]);
 							if (verb)
 								printf("dispcal: inverting Jacobian failed (3) - falling back\n");
 
@@ -4244,6 +4432,7 @@ int main(int argc, char *argv[]) {
 						}
 						/* Restart at the best we've had */
 						if (asgrey.s[i].hde > besthde) {
+							asgrey.s[i]._de = best_de;
 							asgrey.s[i].de = bestde;
 							asgrey.s[i].dc = bestdc;
 							asgrey.s[i].peqde = bestpeqde;
@@ -4257,13 +4446,13 @@ int main(int argc, char *argv[]) {
 							icmSub3(asgrey.s[i].deXYZ, asgrey.s[i].tXYZ, asgrey.s[i].XYZ);
 						}
 						mjac = 1;
-						impj = 1;
+						impj = 1;		/* Have remeasured */
 					}
 #endif /* REMEAS_JACOBIAN */
 
 					/* Compute a correction to the Jacobian if we can. */
 					/* (Don't do this unless we have a solid previous */
-					/* reading for this patch) */ 
+					/* reading for this patch, and we haven't remeasured it) */ 
 					if (impj == 0 && rpt > 0 && isclipped == 0) {
 						double nsdrgb;			/* Norm squared of pdrgb */
 						double spdrgb[3];		/* Scaled previous delta rgb */
@@ -4272,6 +4461,10 @@ int main(int argc, char *argv[]) {
 						double tj[3][3];		/* Temp Jacobian */
 						double itj[3][3];		/* Temp inverse Jacobian */
 
+//printf("~1 Jacobian was: %f %f %f\n", asgrey.s[i].j[0][0], asgrey.s[i].j[0][1], asgrey.s[i].j[0][2]);
+//printf("~1               %f %f %f\n", asgrey.s[i].j[1][0], asgrey.s[i].j[1][1], asgrey.s[i].j[1][2]);
+//printf("~1               %f %f %f\n", asgrey.s[i].j[2][0], asgrey.s[i].j[2][1], asgrey.s[i].j[2][2]);
+
 						/* Use Broyden's Formula */
 						icmSub3(dXYZerr, asgrey.s[i].dXYZ, asgrey.s[i].pdXYZ);
 //printf("~1 Jacobian error = %f %f %f\n", dXYZerr[0], dXYZerr[1], dXYZerr[2]);
@@ -4287,31 +4480,68 @@ int main(int argc, char *argv[]) {
 							{
 								double eXYZ[3];
 	
+//printf("~1 del RGB %f %f %f got del XYZ %f %f %f\n", asgrey.s[i].pdrgb[0], asgrey.s[i].pdrgb[1], asgrey.s[i].pdrgb[2], asgrey.s[i].dXYZ[0], asgrey.s[i].dXYZ[1], asgrey.s[i].dXYZ[2]);
+
 								/* Make a full adjustment to temporary Jac */
 								icmAdd3x3(tj, asgrey.s[i].j, jadj);
+
+//printf("~1 Full Jacobian: %f %f %f\n", tj[0][0], tj[0][1], tj[0][2]);
+//printf("~1                %f %f %f\n", tj[1][0], tj[1][1], tj[1][2]);
+//printf("~1                %f %f %f\n", tj[2][0], tj[2][1], tj[2][2]);
+
 								icmMulBy3x3(eXYZ, tj, asgrey.s[i].pdrgb);
 								icmSub3(eXYZ, eXYZ, asgrey.s[i].dXYZ);
 								printf("Jac check resid %f %f %f\n", eXYZ[0], eXYZ[1], eXYZ[2]);
 							}
 #endif	/* DEBUG */
 
+							/* Add to portion of previous adjustment */
+							/* to counteract undershoot & overshoot */
+							icmScale3x3(pjadj, pjadj, JAC_COMP_FACT);
+							icmAdd3x3(jadj, jadj, pjadj);
+							icmCpy3x3(pjadj, jadj);
+
 							/* Add part of our correction to actual Jacobian */
+							/* to smooth out correction to counteract noise */
 							icmScale3x3(jadj, jadj, JAC_COR_FACT);
 							icmAdd3x3(tj, asgrey.s[i].j, jadj);
+
 							if (icmInverse3x3(itj, tj) == 0) {		/* Invert OK */
 								icmCpy3x3(asgrey.s[i].j, tj);		/* Use adjusted */
 								icmCpy3x3(asgrey.s[i].ij, itj);
 								impj = 1;
+
+#ifdef NEVER
+					/* Check how close new Jacobian predicts previous delta XYZ */
+					{
+						double eXYZ[3];
+						double ergb[3];
+
+						icmMulBy3x3(eXYZ, asgrey.s[i].j, asgrey.s[i].pdrgb);
+						icmSub3(eXYZ, eXYZ, asgrey.s[i].dXYZ);
+						printf("Jac check2 resid %f %f %f\n", eXYZ[0], eXYZ[1], eXYZ[2]);
+
+						icmMulBy3x3(ergb, asgrey.s[i].ij, asgrey.s[i].pdXYZ);
+						printf("Jac check2 drgb would have been %f %f %f\n", ergb[0], ergb[1], ergb[2]);
+						icmAdd3(ergb, ergb, asgrey.s[i].prgb);
+						printf("Jac check2 rgb would have been %f %f %f\n", ergb[0], ergb[1], ergb[2]);
+					}
+#endif
 							}
-//else printf("~1 ij failed\n");
+//else printf("~1 ij failed - reverted\n");
 						}
 //else printf("~1 nsdrgb was below threshold\n");
 					}
 //else if (isclipped) printf("~1 no j update: rgb is clipped\n");
+//printf("~1 Jacobian now: %f %f %f\n", asgrey.s[i].j[0][0], asgrey.s[i].j[0][1], asgrey.s[i].j[0][2]);
+//printf("~1               %f %f %f\n", asgrey.s[i].j[1][0], asgrey.s[i].j[1][1], asgrey.s[i].j[1][2]);
+//printf("~1               %f %f %f\n", asgrey.s[i].j[2][0], asgrey.s[i].j[2][1], asgrey.s[i].j[2][2]);
+
 #endif	/* ADJ_JACOBIAN */
 
 					/* Track the best solution we've found */
 					if (asgrey.s[i].hde <= besthde) {
+						best_de = asgrey.s[i]._de;
 						bestde = asgrey.s[i].de;
 						bestdc = asgrey.s[i].dc;
 						bestpeqde = asgrey.s[i].peqde;
@@ -4323,13 +4553,14 @@ int main(int argc, char *argv[]) {
 						bestxyz[1] = asgrey.s[i].XYZ[1];
 						bestxyz[2] = asgrey.s[i].XYZ[2];
 
+//printf("~1 new best\n");
 					} else if (asgrey.s[i].dc > bestdc) {
 						/* we got worse in Lab and XYZ ! */
 
-						wde = asgrey.s[i].de;
-
 						/* If we've wandered too far, return to best we found */
 						if (asgrey.s[i].hde > (3.0 * besthde)) {
+//printf("~1 resetting to last best\n");
+							asgrey.s[i]._de = best_de;
 							asgrey.s[i].de = bestde;
 							asgrey.s[i].dc = bestdc;
 							asgrey.s[i].peqde = bestpeqde;
@@ -4344,22 +4575,26 @@ int main(int argc, char *argv[]) {
 						}
 
 						/* If the Jacobian hasn't changed, moderate the gain */
-						if (impj == 0)
+						if (impj == 0) {
 							rgain *= 0.8;		/* We might be overshooting */
+//printf("~1 reducing rgain to %f\n",rgain);
+						}
 						gworse = 1;
 					}
 
 					/* See if we need to repeat */
-					if (asgrey.s[i].de <= errthr && asgrey.s[i].peqde < errthr) {	
-						if (verb > 1)
+					if (asgrey.s[i].de <= ierrth && asgrey.s[i].peqde < ierrth) {	
+						if (verb > 1) {
 							if (it < (mxits-1))
-								printf("Point %d Delta E %f, OK\n",rsteps - i,asgrey.s[i].de);
+								printf("Point %d DE %f, W.DE %f, OK ( < %f)\n",rsteps - i,asgrey.s[i]._de, asgrey.s[i].de, ierrth);
 							else
-								printf("Point %d Delta E %f, peqDE %f, OK\n",rsteps - i,asgrey.s[i].de, asgrey.s[i].peqde);
+								printf("Point %d DE %f, W.DE %f, W.peqDE %f, OK ( < %f)\n",rsteps - i,asgrey.s[i]._de,asgrey.s[i].de, asgrey.s[i].peqde, ierrth);
+						}
 						break;	/* No more retries */
 					}
-					if ((rpt+1) >= MAX_RPTS) {
-						asgrey.s[i].de = bestde;			/* Restore to best we found */
+					if ((rpt+1) >= mxrpts) {
+						asgrey.s[i]._de = best_de;			/* Restore to best we found */
+						asgrey.s[i].de = bestde;
 						asgrey.s[i].dc = bestdc;
 						asgrey.s[i].peqde = bestpeqde;		/* Restore to best we found */
 						asgrey.s[i].hde = besthde;			/* Restore to best we found */
@@ -4369,11 +4604,12 @@ int main(int argc, char *argv[]) {
 						asgrey.s[i].XYZ[0] = bestxyz[0];
 						asgrey.s[i].XYZ[1] = bestxyz[1];
 						asgrey.s[i].XYZ[2] = bestxyz[2];
-						if (verb > 1)
+						if (verb > 1) {
 							if (it < (mxits-1))
-								printf("Point %d Delta E %f, Fail\n",rsteps - i,asgrey.s[i].de);
+								printf("Point %d DE %f, W.DE %f, Fail ( > %f)\n",rsteps - i,asgrey.s[i]._de, asgrey.s[i].de, ierrth);
 							else
-								printf("Point %d Delta E %f, peqDE %f, Fail\n",rsteps - i,asgrey.s[i].de,asgrey.s[i].peqde);
+								printf("Point %d DE %f, W.DE %f, W.peqDE %f, Fail ( > %f)\n",rsteps - i,asgrey.s[i]._de,asgrey.s[i].de,asgrey.s[i].peqde,ierrth);
+						}
 						thrfail = 1;						/* Failed to meet target */
 						if (bestde > failerr)
 							failerr = bestde;				/* Worst failed delta E */
@@ -4382,26 +4618,28 @@ int main(int argc, char *argv[]) {
 					if (verb > 1) {
 						if (gworse)
 							if (it < (mxits-1))
-								printf("Point %d Delta E %f, Repeat (got worse)\n", rsteps - i, wde);
+								printf("Point %d DE %f, W.DE %f, Repeat (got worse)\n", rsteps - i, w_de, wde);
 							else
-								printf("Point %d Delta E %f, peqDE %f, Repeat (got worse)\n", rsteps - i, wde,asgrey.s[i].peqde);
+								printf("Point %d DE %f, W.DE %f, peqDE %f, Repeat (got worse)\n", rsteps - i, w_de, wde,asgrey.s[i].peqde);
 						else
 							if (it < (mxits-1))
-								printf("Point %d Delta E %f, Repeat\n", rsteps - i,asgrey.s[i].de);
+								printf("Point %d DE %f, W.DE %f, Repeat\n", rsteps - i,asgrey.s[i]._de,asgrey.s[i].de);
 							else
-								printf("Point %d Delta E %f, peqDE %f, Repeat\n", rsteps - i,asgrey.s[i].de,asgrey.s[i].peqde);
+								printf("Point %d DE %f, W.DE %f, peqDE %f, Repeat\n", rsteps - i,asgrey.s[i]._de,asgrey.s[i].de,asgrey.s[i].peqde);
 					}
 				
+//printf("~1 RGB Jacobian: %f %f %f\n", asgrey.s[i].j[0][0], asgrey.s[i].j[0][1], asgrey.s[i].j[0][2]);
+//printf("~1               %f %f %f\n", asgrey.s[i].j[1][0], asgrey.s[i].j[1][1], asgrey.s[i].j[1][2]);
+//printf("~1               %f %f %f\n", asgrey.s[i].j[2][0], asgrey.s[i].j[2][1], asgrey.s[i].j[2][2]);
 					/* Compute refinement of rgb */
 					icmMulBy3x3(asgrey.s[i].pdrgb, asgrey.s[i].ij, asgrey.s[i].deXYZ);
-//printf("~1 delta needed %f %f %f -> delta RGB %f %f %f\n",
+//printf("~1 XYZ delta needed %f %f %f -> delta RGB %f %f %f\n",
 //asgrey.s[i].deXYZ[0], asgrey.s[i].deXYZ[1], asgrey.s[i].deXYZ[2],
 //asgrey.s[i].pdrgb[0], asgrey.s[i].pdrgb[1], asgrey.s[i].pdrgb[2]);
 
 					/* Gain scale */
 					icmScale3(asgrey.s[i].pdrgb, asgrey.s[i].pdrgb, rgain);
-//printf("~1 delta RGB after gain scale %f %f %f\n",
-//asgrey.s[i].pdrgb[0], asgrey.s[i].pdrgb[1], asgrey.s[i].pdrgb[2]);
+//printf("~1 delta RGB after gain scale %f %f %f\n", asgrey.s[i].pdrgb[0], asgrey.s[i].pdrgb[1], asgrey.s[i].pdrgb[2]);
 
 #ifdef CLIP
 					/* Component wise clip */
@@ -4415,28 +4653,39 @@ int main(int argc, char *argv[]) {
 							dclip = 1;
 						}
 					}
-#ifdef DEBUG
-					if (dclip) printf("delta RGB after clip %f %f %f\n",
+					if (verb >= 3 && dclip) printf("delta RGB after clip %f %f %f\n",
 					       asgrey.s[i].pdrgb[0], asgrey.s[i].pdrgb[1], asgrey.s[i].pdrgb[2]);
-#endif	/* DEBUG */
 #endif	/* CLIP */
 					/* Compute next on the basis of this one RGB */
+					icmCpy3(asgrey.s[i].prgb, asgrey.s[i].rgb);		/* Save previous */
 					icmAdd3(asgrey.s[i].rgb, asgrey.s[i].rgb, asgrey.s[i].pdrgb);
 
 					/* Save expected change in XYZ */
 					icmMulBy3x3(asgrey.s[i].pdXYZ, asgrey.s[i].j, asgrey.s[i].pdrgb);
-#ifdef DEBUG
-					printf("New rgb %f %f %f from expected del XYZ %f %f %f\n",
-					       asgrey.s[i].rgb[0], asgrey.s[i].rgb[1], asgrey.s[i].rgb[2],
-					       asgrey.s[i].pdXYZ[0], asgrey.s[i].pdXYZ[1], asgrey.s[i].pdXYZ[2]);
-#endif
+					if (verb >= 3) {
+						printf("New rgb %f %f %f from expected del XYZ %f %f %f\n",
+						       asgrey.s[i].rgb[0], asgrey.s[i].rgb[1], asgrey.s[i].rgb[2],
+						       asgrey.s[i].pdXYZ[0], asgrey.s[i].pdXYZ[1], asgrey.s[i].pdXYZ[2]);
+					}
 				} else {	/* Verification, so no repeat */
 					break;
 				}
 
 				prevde = asgrey.s[i].de;
 			}	/* Next repeat */
-		}		/* Next resolution step */
+
+			if (verb >= 3) {
+				printf("After adjustment:\n");
+				printf("Current rgb %f %f %f -> XYZ %f %f %f, de %f, dc %f\n", 
+				asgrey.s[i].rgb[0], asgrey.s[i].rgb[1], asgrey.s[i].rgb[2],
+				asgrey.s[i].XYZ[0], asgrey.s[i].XYZ[1], asgrey.s[i].XYZ[2],
+				asgrey.s[i].de, asgrey.s[i].dc);
+				printf("Target XYZ %f %f %f, delta needed %f %f %f\n", 
+				asgrey.s[i].tXYZ[0], asgrey.s[i].tXYZ[1], asgrey.s[i].tXYZ[2],
+				asgrey.s[i].deXYZ[0], asgrey.s[i].deXYZ[1], asgrey.s[i].deXYZ[2]);
+			}
+
+		}		/* Next patch/step */
 		if (verb)
 			printf("\n");			/* Final return for patch count */
 
@@ -4487,7 +4736,7 @@ int main(int argc, char *argv[]) {
 			/* We're checking against our given brightness and */
 			/* white point target. */
 			mnerr = anerr = 0.0;
-			init_csamp_txyz(&asgrey, &x, 0);		/* In case the targets were tweaked */
+			init_csamp_txyz(&asgrey, &x, 0, verb);	/* In case the targets were tweaked */
 			for (i = 0; i < asgrey.no; i++) {
 				double err;
 
@@ -4518,7 +4767,7 @@ int main(int argc, char *argv[]) {
 		}
 
 		/* Verify loop exit */
-		if (it >= (mxits + nver -1)) {
+		if (nver > 0 && it >= (mxits + nver -1)) {
 			break;
 		}
 
@@ -4546,7 +4795,6 @@ int main(int argc, char *argv[]) {
 					sdv[j][i].p = asgrey.s[i].v;
 					sdv[j][i].v = asgrey.s[i].rgb[j];
 					sdv[j][i].w = 1.0;
-// ~~999
 #ifdef NEVER
 					printf("rdac %d point %d = %f, %f\n",j,i,sdv[j][i].p,sdv[j][i].v);
 #endif
@@ -4623,6 +4871,8 @@ int main(int argc, char *argv[]) {
 		/* Tell downstream whether they can expect that this calibration */
 		/* will be applied in hardware or not. */
 		ocg->add_kword(ocg, 0, "VIDEO_LUT_CALIBRATION_POSSIBLE",noramdac ? "NO" : "YES", NULL);
+		/* Tell downstream whether the device range was actually (16-235)/255 */
+		ocg->add_kword(ocg, 0, "TV_OUTPUT_ENCODING",out_tvenc ? "YES" : "NO", NULL);
 
 		/* Put the target parameters in the CGATS file too */
 		if (dtype != 0) {
@@ -4833,6 +5083,16 @@ int main(int argc, char *argv[]) {
 					cc = 0.0;
 				else if (cc > 1.0)
 					cc = 1.0;
+				if (out_tvenc) {
+					cc = (cc * (235.0-16.0) + 16.0)/255.0;
+
+					/* For video encoding the extra bits of precision are created by bit shifting */
+					/* rather than scaling, so we need to scale the fp value to account for this. */
+					/* We assume the precision is the vcgt table size = 16 */
+					/* ~~99 ideally we should tag the fact that this is video encoded, so that */
+					/* the vcgt loaded can adjust for a different bit precision ~~~~ */
+					cc = (cc * 255 * (1 << (16 - 8)))/((1 << 16) - 1.0); 	
+				}
 				((unsigned short*)wo->u.table.data)[CAL_RES * j + i] = (int)(cc * 65535.0 + 0.5);
 			}
 		}
@@ -5174,6 +5434,15 @@ int main(int argc, char *argv[]) {
 						cc = 0.0;
 					else if (cc > 1.0)
 						cc = 1.0;
+					if (out_tvenc) {
+						cc = (cc * (235.0-16.0) + 16.0)/255.0;
+						/* For video encoding the extra bits of precision are created by bit */
+						/* shifting rather than scaling, so we need to scale the fp value to */
+						/* account for this. We assume the precision is the vcgt table size = 16 */
+						/* ~~99 ideally we should tag the fact that this is video encoded, so */
+						/* that the vcgt loaded can adjust for a different bit precision ~~~~ */
+						cc = (cc * 255 * (1 << (16 - 8)))/((1 << 16) - 1.0); 	
+					}
 					((unsigned short*)wo->u.table.data)[CAL_RES * j + i] = (int)(cc * 65535.0 + 0.5);
 				}
 			}
@@ -5197,9 +5466,13 @@ int main(int argc, char *argv[]) {
 			wog->allocate((icmBase *)wog);
 			wob->allocate((icmBase *)wob);
 
-			wor->data[0].X = mat[0][0]; wor->data[0].Y = mat[1][0]; wor->data[0].Z = mat[2][0];
-			wog->data[0].X = mat[0][1]; wog->data[0].Y = mat[1][1]; wog->data[0].Z = mat[2][1];
-			wob->data[0].X = mat[0][2]; wob->data[0].Y = mat[1][2]; wob->data[0].Z = mat[2][2];
+			/* Make sure rounding doesn't wreck white point */
+			icmTranspose3x3(mat, mat);	/* Convert [XYZ][RGB] to [RGB][XYZ] */
+			quantizeRGBprimsS15Fixed16(mat);
+
+			wor->data[0].X = mat[0][0]; wor->data[0].Y = mat[0][1]; wor->data[0].Z = mat[0][2];
+			wog->data[0].X = mat[1][0]; wog->data[0].Y = mat[1][1]; wog->data[0].Z = mat[1][2];
+			wob->data[0].X = mat[2][0]; wob->data[0].Y = mat[2][1]; wob->data[0].Z = mat[2][2];
 		}
 
 		/* Red, Green and Blue Tone Reproduction Curve Tags: */