Imported Upstream version 1.7.0upstream/1.7.0

45 files changed, 1561 insertions, 1250 deletions

diff --git a/xicc/Jamfile b/xicc/Jamfile
index 58e5051..9c15349 100644
--- a/xicc/Jamfile
+++ b/xicc/Jamfile
@@ -21,7 +21,7 @@ HDRS = ../h ../icc ../rspl ../cgats ../numlib ../gamut ../spectro ../profile
 
 # XICC library
 Library libxicc :  xicc.c xlutfix.c xspect.c xcolorants.c xutils.c iccjpeg.c xdevlin.c
-        xcam.c cam97s3.c cam02.c mpp.c ccmx.c ccss.c xfit.c xdgb.c moncurve.c xcal.c ;
+        xcam.c cam97s3.c cam02.c mpp.c ccmx.c ccss.c xfit.c xdgb.c moncurve.c xcal.c bt1886.c ;
 
 # colorant library. Use instead of libxicc
 Object xcolorants2 : xcolorants.c ;
@@ -35,7 +35,8 @@ LibraryFromObjects libxutils : xutils2 iccjpeg2 ;
 # Utilities / test programs
 
 LINKLIBS = libxicc ../spectro/libinsttypes ../gamut/libgamut ../rspl/librspl
-           ../cgats/libcgats ../icc/libicc ../plot/libplot ../plot/libvrml ../numlib/libnum
+           ../cgats/libcgats ../icc/libicc ../plot/libplot ../plot/libvrml
+           ../numlib/libnum ../numlib/libui
            $(TIFFLIB) $(JPEGLIB) ; 
 
 # Not created yet
@@ -54,7 +55,7 @@ Main iccgamut : iccgamut.c ;
 Main tiffgamut : tiffgamut.c : : : $(TIFFINC) $(JPEGINC) : : ;
 
 # diagnostic utility
-if [ GLOB . : tiffgmts.c ]  {
+if [ GLOB [ NormPaths . ] : tiffgmts.c ]  {
 	Main tiffgmts : tiffgmts.c : : : $(TIFFINC) $(JPEGINC) : : ../plot/libvrml ;
 }
 
@@ -137,6 +138,8 @@ if $(BUILD_JUNK) {
 
 	LINKLIBS += ../render/librender ;
 
+	Main illlocus : illlocus.c ;
+
 	Main slocustest : slocustest.c ;
 	
 	MainsFromSources t1.c t2.c t22.c t23.c t24.c t3.c ;
diff --git a/xicc/afiles b/xicc/afiles
index 8b4341a..f2b84ad 100644
--- a/xicc/afiles
+++ b/xicc/afiles
@@ -54,6 +54,8 @@ cam02plot.c
 moncurve.c
 moncurve.h
 monctest.c
+bt1886.h
+bt1886.c
 fbview.c
 xfbview.c
 icheck.c
diff --git a/xicc/bt1886.c b/xicc/bt1886.c
new file mode 100644
index 0000000..3f94dbe
--- /dev/null
+++ b/xicc/bt1886.c
@@ -0,0 +1,351 @@
+
+/* 
+ * Author:  Graeme W. Gill
+ * Date:    16/8/13
+ * Version: 1.00
+ *
+ * Copyright 2013, 2014 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU AFFERO GENERAL PUB LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ *
+ */
+
+/* BT.1886 stype input offset transfer curve, */
+/* + general gamma + input + output offset curve support. */
+
+#include <sys/types.h>
+#include <string.h>
+#include <ctype.h>
+#ifdef __sun
+#include <unistd.h>
+#endif
+#if defined(__IBMC__) && defined(_M_IX86)
+#include <float.h>
+#endif
+#include "numlib.h"
+#include "icc.h"		/* definitions for this library */
+#include "bt1886.h"		/* definitions for this library */
+
+#undef DEBUG
+
+/* BT.1886 support */
+/* This is both a EOTF curve, and a white point */
+/* adjustment. */
+
+/* Compute technical gamma from effective gamma in BT.1886 style */
+
+/* Info for optimization */
+typedef struct {
+	double wp;			/* 100% input target */ 
+	double thyr;		/* 50% input target */
+	double bp;			/* 0% input target */
+} gam_fits;
+
+/* gamma + input offset function handed to powell() */
+static double gam_fit(void *dd, double *v) {
+	gam_fits *gf = (gam_fits *)dd;
+	double gamma = v[0];
+	double a, b;
+	double rv = 0.0;
+	double t1, t2;
+
+	if (gamma < 0.0) {
+		rv += 100.0 * -gamma;
+		gamma = 1e-4;
+	}
+
+	t1 = pow(gf->bp, 1.0/gamma);
+	t2 = pow(gf->wp, 1.0/gamma);
+	b = t1/(t2 - t1);						/* Offset */
+	a = pow(t2 - t1, gamma);				/* Gain */
+
+	/* Comput 50% output for this technical gamma */
+	/* (All values are without output offset being added in) */
+	t1 = a * pow((0.5 + b), gamma);
+	t1 = t1 - gf->thyr;
+	rv += t1 * t1;
+	
+	return rv;
+}
+
+/* Given the effective gamma and the output offset Y, */
+/* return the technical gamma needed for the correct 50% response. */
+static double xicc_tech_gamma(
+	double egamma,			/* effective gamma needed */
+	double off,				/* Output offset required */
+	double outoprop			/* Prop. of offset to be accounted for on output */
+) {
+	gam_fits gf;
+	double outo;
+	double op[1], sa[1], rv;
+
+	if (off <= 0.0) {
+		return egamma;
+	}
+
+	/* We set up targets without outo being added */
+	outo = off * outoprop; 		/* Offset acounted for in output */
+	gf.bp = off - outo;			/* Black value for 0 % input */
+	gf.wp = 1.0 - outo;			/* White value for 100% input */
+	gf.thyr = pow(0.5, egamma) - outo;	/* Advetised 50% target */
+
+	op[0] = egamma;
+	sa[0] = 0.1;
+
+	if (powell(&rv, 1, op, sa, 1e-6, 500, gam_fit, (void *)&gf, NULL, NULL) != 0)
+		warning("Computing effective gamma and input offset is inaccurate");
+
+	return op[0];
+}
+
+/* Set the bt1886_info to a default do nothing state */
+void bt1886_setnop(bt1886_info *p) {
+	icmXYZ2XYZ(p->w, icmD50);
+	p->ingo = 0.0;
+	p->outsc = 1.0;
+	p->outo = 0.0;
+	p->outL = 0.0;
+	p->tab[1] = 0.0;
+	p->tab[2] = 0.0;
+}
+
+/* Setup the bt1886_info for the given target black point, proportion of */
+/* offset to be accounted for on output, and gamma. */
+/* wp XYZ simply sets the L*a*b* reference */
+/* Pure BT.1886 will have outopro = 0.0 and gamma = 2.4 */
+void bt1886_setup(
+bt1886_info *p,
+icmXYZNumber *w,	/* wp used for L*a*b* conversion */
+double *XYZbp,		/* normalised bp used for black offset and black point hue "bend" */
+double outoprop,	/* 0..1 proportion of output black point compensation */
+double gamma,		/* technical or effective gamma */
+int effg			/* nz if effective gamma, z if technical gamma */
+) {
+	double Lab[3], ino, bkipow, wtipow;
+
+	icmXYZ2XYZ(p->w, *w);
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886_setup wp.Y %f, bp.Y %f, outprop %f, gamma %f, effg %d", p->w.Y, XYZbp[1], outoprop, gamma, effg);
+#endif
+
+	if (effg) {
+		p->gamma = xicc_tech_gamma(gamma, XYZbp[1], outoprop);
+#ifdef DEBUG
+		a1logd(g_log, 2, "bt1886_setup tgamma %f", p->gamma);
+#endif
+	} else {
+		p->gamma = gamma;
+	}
+
+	icmXYZ2Lab(&p->w, Lab, XYZbp);
+
+	p->outL   = Lab[0];		/* For bp blend comp. */
+	p->tab[0] = Lab[1];		/* a* b* correction needed */
+	p->tab[1] = Lab[2];
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886_setup bend Lab = %f %f %f", p->outL, p->tab[0], p->tab[1]);
+#endif
+
+	if (XYZbp[1] < 0.0)
+		XYZbp[1] = 0.0;
+
+
+	p->outo = XYZbp[1] * outoprop; 		/* Offset acounted for in output */
+	ino = XYZbp[1] - p->outo;			/* Balance of offset accounted for in input */
+
+	bkipow = pow(ino, 1.0/p->gamma);				/* Input offset black to 1/pow */
+	wtipow = pow((1.0 - p->outo), 1.0/p->gamma);	/* Input offset white to 1/pow */
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886_setup outo %f, ino %f, bkipow %f, wtipow %f", p->outo, ino, bkipow, wtipow);
+#endif
+
+	p->ingo = bkipow/(wtipow - bkipow);			/* non-linear Y that makes input offset */
+												/* proportion of black point */
+	p->outsc = pow(wtipow - bkipow, p->gamma);	/* Scale to make input of 1 map to */
+												/* 1.0 - p->outo */
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886_setup ingo %f, outsc %f", p->ingo, p->outsc);
+#endif
+}
+
+/* Apply BT.1886 eotf curve to the device RGB value */
+/* to produce a linear light RGB. We pass xvYCC out of range values through. */
+void bt1886_fwd_curve(bt1886_info *p, double *out, double *in) {
+	int j;
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886    Dev RGB in %f %f %f, pow %f\n", in[0],in[1],in[2], p->gamma);
+	a1logd(g_log, 2, "outo %f, outsc %f, pow %f\n", p->outo, p->outsc, 1.0/p->gamma);
+	a1logd(g_log, 2, "ingo %f, pow %f, outsc %f, outo %f\n", p->ingo, p->gamma, p->outsc,p->outo);
+#endif
+
+	for (j = 0; j < 3; j++) {
+		int neg = 0;
+		double vv = in[j];
+	
+		if (vv < 0.0) {			/* Allow for xvYCC */
+			neg = 1;
+			vv = -vv;
+		}
+		/* Apply input offset */
+		vv += p->ingo;
+		
+		/* Apply power and scale */
+		if (vv > 0.0)
+			vv = p->outsc * pow(vv, p->gamma);
+
+		/* Apply output portion of offset */
+		vv += p->outo;
+
+		if (neg)
+			vv = -vv;
+
+		out[j] = vv;
+	}
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 linear RGB out %f %f %f\n", out[0],out[1],out[2]);
+#endif
+}
+
+/* Apply inverse BT.1886 eotf curve to the linear light RGB to produce */
+/* device RGB values. We pass xvYCC out of range values through. */
+void bt1886_bwd_curve(bt1886_info *p, double *out, double *in) {
+	int j;
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 linear RGB in %f %f %f\n", in[0],in[1],in[2]);
+	a1logd(g_log, 2, "outo %f, outsc %f, pow %f, ingo %f\n", p->outo, p->outsc, 1.0/p->gamma,p->ingo);
+#endif
+
+	for (j = 0; j < 3; j++) {
+		int neg = 0;
+		double vv = in[j];
+	
+		if (vv < 0.0) {			/* Allow for xvYCC */
+			neg = 1;
+			vv = -vv;
+		}
+
+		/* Un-apply output portion of offset */
+		vv -= p->outo;
+
+		/* Un-apply power and scale */
+		if (vv > 0.0)
+			vv = pow(vv/p->outsc, 1.0/p->gamma);
+
+		/* Un-apply input offset */
+		vv -= p->ingo;
+		
+		if (neg)
+			vv = -vv;
+
+		out[j] = vv;
+	}
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886   Dev RGB out %f %f %f\n", in[0],in[1],in[2]);
+#endif
+}
+
+/* Apply BT.1886 processing black point hue adjustment to the XYZ value */
+void bt1886_wp_adjust(bt1886_info *p, double *out, double *in) {
+	double vv;
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 XYZ wp adj. in %f %f %f\n", in[0],in[1],in[2]);
+#endif
+
+	icmXYZ2Lab(&p->w, out, in);
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 Lab wp adj. in %f %f %f\n", out[0],out[1],out[2]);
+#endif
+
+	/* Blend ab to required black point offset p->tab[] as L approaches black. */
+	vv = (out[0] - p->outL)/(100.0 - p->outL);	/* 0 at bp, 1 at wp */
+	vv = 1.0 - vv;
+	if (vv < 0.0)
+		vv = 0.0;
+	else if (vv > 1.0)
+		vv = 1.0;
+	vv = pow(vv, 40.0);
+
+	out[1] += vv * p->tab[0];
+	out[2] += vv * p->tab[1];
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 Lab after wp adj. %f %f %f\n", out[0],out[1],out[2]);
+#endif
+
+	icmLab2XYZ(&p->w, out, out);
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 XYZ after wp adj. %f %f %f\n", out[0],out[1],out[2]);
+#endif
+}
+
+
+/* Apply inverse BT.1886 processing black point hue adjustment to the XYZ value */
+void bt1886_inv_wp_adjust(bt1886_info *p, double *out, double *in) {
+	double vv;
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 XYZ inv. wp adj. in %f %f %f\n", in[0],in[1],in[2]);
+#endif
+
+	icmXYZ2Lab(&p->w, out, in);
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 Lab inv. wp adj. in %f %f %f\n", out[0],out[1],out[2]);
+#endif
+
+	/* Blend ab to required black point offset p->tab[] as L approaches black. */
+	vv = (out[0] - p->outL)/(100.0 - p->outL);	/* 0 at bp, 1 at wp */
+	vv = 1.0 - vv;
+	if (vv < 0.0)
+		vv = 0.0;
+	else if (vv > 1.0)
+		vv = 1.0;
+	vv = pow(vv, 40.0);
+
+	out[1] -= vv * p->tab[0];
+	out[2] -= vv * p->tab[1];
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 Lab after inv. wp adj. %f %f %f\n", out[0],out[1],out[2]);
+#endif
+
+	icmLab2XYZ(&p->w, out, out);
+
+#ifdef DEBUG
+	a1logd(g_log, 2, "bt1886 XYZ after inv. wp adj. %f %f %f\n", out[0],out[1],out[2]);
+#endif
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/xicc/bt1886.h b/xicc/bt1886.h
new file mode 100644
index 0000000..f14efdc
--- /dev/null
+++ b/xicc/bt1886.h
@@ -0,0 +1,94 @@
+#ifndef BT1886_H
+#define BT1886_H
+
+/* 
+ * Author:  Graeme W. Gill
+ * Date:    16/8/13
+ * Version: 1.00
+ *
+ * Copyright 2013, 2014 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU AFFERO GENERAL PUB LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ *
+ */
+
+/* BT.1886 stype input offset transfer curve, */
+/* + general gamma + input + output offset curve support. */
+
+
+typedef struct {
+	icmXYZNumber w;				/* White point for Lab conversion */
+	double ingo;				/* input Y offset for bt1886 */
+	double outsc;				/* output Y scale for bt1886 */
+	double outo;				/* output Y offset */
+	double outL;				/* output black point L value */
+	double tab[2];				/* Target ab offset value at zero input for bt1886 */
+	double gamma;				/* bt.1886 technical gamma to apply */
+} bt1886_info;
+
+/* Set the bt1886_info to a default do nothing state */
+void bt1886_setnop(bt1886_info *p);
+
+/* Setup the bt1886_info for the given target */
+void bt1886_setup(
+bt1886_info *p,
+icmXYZNumber *w,	/* wp used for L*a*b* conversion */
+double *XYZbp,		/* normalised bp used for black offset and black point hue "bend" */
+double outoprop,	/* 0..1 proportion of output black point compensation */
+double gamma,		/* technical or effective gamma */
+int effg			/* nz if effective gamma, z if technical gamma */
+);
+
+/* Apply BT.1886 eotf curve to the device RGB value */
+/* to produce a linear light RGB. We pass xvYCC out of range values through. */
+void bt1886_fwd_curve(bt1886_info *p, double *out, double *in);
+
+/* Apply BT.1886 processing black point hue adjustment to the XYZ value */
+void bt1886_wp_adjust(bt1886_info *p, double *out, double *in);
+
+
+/* Apply inverse BT.1886 eotf curve to the device RGB value */
+/* to produce a linear light RGB. We pass xvYCC out of range values through. */
+void bt1886_bwd_curve(bt1886_info *p, double *out, double *in);
+
+/* Apply inverse BT.1886 processing black point hue adjustment to the XYZ value */
+void bt1886_inv_wp_adjust(bt1886_info *p, double *out, double *in);
+
+#endif /* BT1886_H */
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/xicc/cam02.c b/xicc/cam02.c
index 3684a78..9bc9343 100644
--- a/xicc/cam02.c
+++ b/xicc/cam02.c
@@ -46,7 +46,7 @@
 	uniformity. A color with one component near zero might shift
 	all the components to -ve values on inverse conversion - ie.
 	a 1 DE shift in Jab becomes a masive DE in XYZ/Lab/perceptual,
-	with (say) a darl red becomong black because the blue
+	with (say) a dark red becomong black because the blue
 	value is small. One way around this is to re-introduce
 	a flag to turn off perfect symetry by disabling
 	expansion on the reverse conversion.
@@ -279,7 +279,7 @@ double Yb,		/* Relative Luminance of Background to reference white (range 0.0 ..
 double Lv,		/* Luminance of white in the Viewing/Scene/Image field (cd/m^2) */
 				/* Ignored if Ev is set to other than vc_none */
 double Yf,		/* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
-double Yg,		/* Flare as a fraction of the ambient (Y range 0.0 .. 1.0) */
+double Yg,		/* Flare as a fraction of the adapting/surround (Y range 0.0 .. 1.0) */
 double Gxyz[3],	/* The Glare white coordinates (typically the Ambient color) */
 				/* If <= 0 will Wxyz will be used. */
 int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
@@ -403,7 +403,7 @@ int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
 	s->Fsxyz[1] = s->Yf * s->Wxyz[1];
 	s->Fsxyz[2] = s->Yf * s->Wxyz[2];
 
-	/* Add in the Glare contribution from the ambient */
+	/* Add in the Glare contribution from the adapting/surround */
 	tt = s->Yg * s->La/s->Lv;
 	s->Fsxyz[0] += tt * s->Gxyz[0];
 	s->Fsxyz[1] += tt * s->Gxyz[1];
@@ -558,8 +558,8 @@ int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
 	printf("Scene parameters:\n");
 	printf("Viewing condition Ev = %d\n",s->Ev);
 	printf("Ref white Wxyz = %f %f %f\n", s->Wxyz[0], s->Wxyz[1], s->Wxyz[2]);
-	printf("Relative liminance of background Yb = %f\n", s->Yb);
-	printf("Adapting liminance La = %f\n", s->La);
+	printf("Relative luminance of background Yb = %f\n", s->Yb);
+	printf("Adapting luminance La = %f\n", s->La);
 	printf("Flare Yf = %f\n", s->Yf);
 	printf("Glare Yg = %f\n", s->Yg);
 	printf("Glare color Gxyz = %f %f %f\n", s->Gxyz[0], s->Gxyz[1], s->Gxyz[2]);
diff --git a/xicc/cam02.h b/xicc/cam02.h
index a42e71f..37443d3 100644
--- a/xicc/cam02.h
+++ b/xicc/cam02.h
@@ -17,7 +17,7 @@
  *
  * This file is based on cam97s3.h by Graeme Gill.
  *
- * Copyright 2004 - 2013 Graeme W. Gill
+ * Copyright 2004 - 2014 Graeme W. Gill
  * Please refer to COPYRIGHT file for details.
  * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
  * see the License.txt file for licencing details.
@@ -39,7 +39,8 @@
   Surround/Adapting field is the visual field minus the background field,
   and is what is assumed to be setting the viewers light adaptation level.
 
-  Ambient field is the whole surrounding environmental light field.
+  Ambient field is general illumination from the sun or general lighting,
+  and is the whole surrounding environmental light field.
 
   Illuminating field is the field that illuminates the reflective
   Scene/Image. It may be the same as the Ambient field or it could
@@ -88,7 +89,7 @@
 
 /* Flare is assumed to be stray light from light parts of the */
 /* image illuminating dark parts of the image, and is display technology dependent. */
-/* In theory reflective systems have no Flare ? */
+/* In theory reflective systems have little Flare ? */
 
 /* Glare is assumed to be stray ambient light reflecting from the display */
 /* surface, dust, or entering the observers eye directly, and as a result */
@@ -138,7 +139,7 @@ struct _cam02 {
 		double Lv,		/* Luminance of white in the Viewing/Scene/Image field (cd/m^2) */
 						/* Ignored if Ev is set */
 		double Yf,		/* Flare as a fraction of the reference white (range 0.0 .. 1.0) */
-		double Yg,		/* Glare as a fraction of the ambient (range 0.0 .. 1.0) */
+		double Yg,		/* Glare as a fraction of the adapting/surround (range 0.0 .. 1.0) */
 		double Gxyz[3],	/* The Glare white coordinates (ie. the Ambient color) */
 						/* If <= 0 will Wxyz will be used. */
 		int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
@@ -156,7 +157,7 @@ struct _cam02 {
 	double Wxyz[3];	/* Reference/Adapted White XYZ (Y range 0.0 .. 1.0) */
 	double Yb;		/* Relative Luminance of Background to reference white (Y range 0.0 .. 1.0) */
 	double Yf;		/* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
-	double Yg;		/* Glare as a fraction of the ambient (Y range 0.0 .. 1.0) */
+	double Yg;		/* Glare as a fraction of the adapting/surround (Y range 0.0 .. 1.0) */
 	double Gxyz[3];	/* The Glare white coordinates (typically the Ambient color) */
 
 	/* Internal parameters */
diff --git a/xicc/cam02ref.h b/xicc/cam02ref.h
index 8965b7a..29f7d2a 100644
--- a/xicc/cam02ref.h
+++ b/xicc/cam02ref.h
@@ -41,7 +41,7 @@ struct _cam02ref {
 		double Lv,		/* Luminance of white in the Viewing/Scene/Image field (cd/m^2) */
 						/* Ignored if Ev is set */
 		double Yf,		/* Flare as a fraction of the reference white (range 0.0 .. 1.0) */
-		double Yg,		/* Glare as a fraction of the ambient (range 0.0 .. 1.0) */
+		double Yg,		/* Glare as a fraction of the adapting/surround (range 0.0 .. 1.0) */
 		double Gxyz[3],	/* The Glare white coordinates (typically the Ambient color) */
 		int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
 	);
@@ -158,7 +158,7 @@ double Yb,		/* Relative Luminence of Background to reference white */
 double Lv,		/* Luminence of white in the Viewing/Scene/Image field (cd/m^2) */
 				/* Ignored if Ev is set to other than vc_none */
 double Yf,		/* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
-double Yg,		/* Glare as a fraction of the ambient (Y range 0.0 .. 1.0) */
+double Yg,		/* Glare as a fraction of the adapting/surround (Y range 0.0 .. 1.0) */
 double Gxyz[3],	/* The Glare white coordinates (typically the Ambient color) */
 int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
 ) {
@@ -223,7 +223,7 @@ int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
 	s->Fsxyz[1] = s->Yf * s->Wxyz[1];
 	s->Fsxyz[2] = s->Yf * s->Wxyz[2];
 
-	/* Add in the Glare contribution from the ambient */
+	/* Add in the Glare contribution from the adapting/surround */
 	tt = s->Yg * s->La/s->Lv;
 	s->Fsxyz[0] += tt * s->Gxyz[0];
 	s->Fsxyz[1] += tt * s->Gxyz[1];
diff --git a/xicc/ccmx.c b/xicc/ccmx.c
index b1244af..464b89c 100644
--- a/xicc/ccmx.c
+++ b/xicc/ccmx.c
@@ -37,6 +37,7 @@
 #include "conv.h"
 #endif
 #include "cgats.h"
+#include "disptechs.h"
 #include "ccmx.h"
 
 #ifdef NT       /* You'd think there might be some standards.... */
@@ -78,18 +79,11 @@ cgats **pocg        /* return CGATS structure */
 	ocg->add_kword(ocg, 0, "INSTRUMENT",p->inst, NULL);
 	if (p->disp != NULL)
 		ocg->add_kword(ocg, 0, "DISPLAY",p->disp, NULL);
-	if (p->tech != NULL)
-		ocg->add_kword(ocg, 0, "TECHNOLOGY",p->tech, NULL);
-	if (p->disp == NULL && p->tech == NULL) {
-#ifdef DEBUG
-		fprintf(stdout, "write_ccmx: ccmx doesn't contain display or techology strings");
-#endif
-		sprintf(p->err, "write_ccmx: ccmx doesn't contain display or techology strings");
-		ocg->del(ocg);
-		return 1;
-	}
-	if (p->cbid != 0) {
-		sprintf(buf, "%d", p->cbid);
+
+	ocg->add_kword(ocg, 0, "TECHNOLOGY", disptech_get_id(p->dtech)->strid,NULL);
+
+	if (p->cc_cbid != 0) {
+		sprintf(buf, "%d", p->cc_cbid);
 		ocg->add_kword(ocg, 0, "DISPLAY_TYPE_BASE_ID", buf, NULL);
 	}
 	if (p->refrmode >= 0)
@@ -199,62 +193,54 @@ cgats *icg			/* input cgats structure */
 
 	if (icg->ntables == 0 || icg->t[0].tt != tt_other || icg->t[0].oi != 0) {
 		sprintf(p->err, "read_ccmx: Input file isn't a CCMX format file");
-		icg->del(icg);
 		return 1;
 	}
 	if (icg->ntables != 1) {
 		sprintf(p->err, "Input file doesn't contain exactly one table");
-		icg->del(icg);
 		return 1;
 	}
 	if ((ti = icg->find_kword(icg, 0, "COLOR_REP")) < 0) {
 		sprintf(p->err, "read_ccmx: Input file doesn't contain keyword COLOR_REP");
-		icg->del(icg);
 		return 1;
 	}
 
 	if (strcmp(icg->t[0].kdata[ti],"XYZ") != 0) { 
 		sprintf(p->err, "read_ccmx: Input file doesn't have COLOR_REP of XYZ");
-		icg->del(icg);
 		return 1;
 	}
 
 	if ((ti = icg->find_kword(icg, 0, "DESCRIPTOR")) >= 0) {
 		if ((p->desc = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccmx: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
 
 	if ((ti = icg->find_kword(icg, 0, "INSTRUMENT")) < 0) {
 		sprintf(p->err, "read_ccmx: Input file doesn't contain keyword INSTRUMENT");
-		icg->del(icg);
 		return 1;
 	}
 	if ((p->inst = strdup(icg->t[0].kdata[ti])) == NULL) {
 		sprintf(p->err, "read_ccmx: malloc failed");
-		icg->del(icg);
 		return 2;
 	}
 
 	if ((ti = icg->find_kword(icg, 0, "DISPLAY")) >= 0) {
 		if ((p->disp = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccmx: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
 	if ((ti = icg->find_kword(icg, 0, "TECHNOLOGY")) >= 0) {
 		if ((p->tech = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccmx: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
+		/* Get disptech enum from standard TECHNOLOGY string */
+		p->dtech = disptech_get_strid(p->tech)->dtech;
 	}
 	if (p->disp == NULL && p->tech == NULL) {
 		sprintf(p->err, "read_ccmx: Input file doesn't contain keyword DISPLAY or TECHNOLOGY");
-		icg->del(icg);
 		return 1;
 	}
 	if ((ti = icg->find_kword(icg, 0, "DISPLAY_TYPE_REFRESH")) >= 0) {
@@ -266,15 +252,14 @@ cgats *icg			/* input cgats structure */
 		p->refrmode = -1;
 	}
 	if ((ti = icg->find_kword(icg, 0, "DISPLAY_TYPE_BASE_ID")) >= 0) {
-		p->cbid = atoi(icg->t[0].kdata[ti]);
+		p->cc_cbid = atoi(icg->t[0].kdata[ti]);
 	} else {
-		p->cbid = 0;
+		p->cc_cbid = 0;
 	}
 
 	if ((ti = icg->find_kword(icg, 0, "UI_SELECTORS")) >= 0) {
 		if ((p->sel = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccmx: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
@@ -282,7 +267,6 @@ cgats *icg			/* input cgats structure */
 	if ((ti = icg->find_kword(icg, 0, "REFERENCE")) >= 0) {
 		if ((p->ref = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccmx: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
@@ -291,19 +275,16 @@ cgats *icg			/* input cgats structure */
 	for (i = 0; i < 3; i++) {	/* XYZ fields */
 		if ((spi[i] = icg->find_field(icg, 0, xyzfname[i])) < 0) {
 			sprintf(p->err, "read_ccmx: Input file doesn't contain field %s", xyzfname[i]);
-			icg->del(icg);
 			return 1;
 		}
 		if (icg->t[0].ftype[spi[i]] != r_t) {
 			sprintf(p->err, "read_ccmx: Input file field %s is wrong type", xyzfname[i]);
-			icg->del(icg);
 			return 1;
 		}
 	}
 
 	if (icg->t[0].nsets != 3) {
 		sprintf(p->err, "read_ccmx: Input file doesn't have exactly 3 sets");
-		icg->del(icg);
 		return 1;
 	}
 
@@ -405,8 +386,8 @@ double *in					/* Input XYZ */
 static int set_ccmx(ccmx *p,
 char *desc,			/* General description (optional) */
 char *inst,			/* Instrument description to copy from */
-char *disp,			/* Display make and model (optional if tech) */
-char *tech,			/* Display technology description (optional if disp) */
+char *disp,			/* Display make and model (optional) */
+disptech dtech,		/* Display technology enum */
 int refrmode,		/* Display refresh mode, -1 = unknown, 0 = n, 1 = yes */
 int cbid,			/* Display type calibration base ID, 0 = unknown */
 char *sel,			/* UI selector characters - NULL for none */
@@ -425,13 +406,9 @@ double mtx[3][3]	/* Transform matrix to copy from */
 		sprintf(p->err, "set_ccmx: malloc failed");
 		return 2;
 	}
-	if ((p->tech = tech) != NULL && (p->tech = strdup(tech)) == NULL) {
-		sprintf(p->err, "set_ccmx: malloc failed");
-		return 2;
-	}
-
+	p->dtech = dtech;
 	p->refrmode = refrmode;
-	p->cbid = cbid;
+	p->cc_cbid = cbid;
 
 	if (sel != NULL) {
 		if ((p->sel = strdup(sel)) == NULL) {
@@ -557,8 +534,8 @@ double optf(void *fdata, double *tp) {
 static int create_ccmx(ccmx *p,
 char *desc,			/* General description (optional) */
 char *inst,			/* Instrument description to copy from */
-char *disp,			/* Display make and model (optional if tech) */
-char *tech,			/* Display technology description (optional if disp) */
+char *disp,			/* Display make and model (optional) */
+disptech dtech,		/* Display technology enum */
 int refrmode,		/* Display refresh mode, -1 = unknown, 0 = n, 1 = yes */
 int cbid,			/* Display type calibration base index, 0 = unknown */
 char *sel,			/* UI selector characters - NULL for none */
@@ -584,12 +561,9 @@ double (*cols)[3]		/* Array of XYZ values from colorimeter */
 		sprintf(p->err, "create_ccmx: malloc failed");
 		return 2;
 	}
-	if ((p->tech = tech) != NULL && (p->tech = strdup(tech)) == NULL) {
-		sprintf(p->err, "create_ccmx: malloc failed");
-		return 2;
-	}
+	p->dtech = dtech;
 	p->refrmode = refrmode;
-	p->cbid = cbid;
+	p->cc_cbid = cbid;
 
 	if (sel != NULL) {
 		if ((p->sel = strdup(sel)) == NULL) {
@@ -680,10 +654,10 @@ double (*cols)[3]		/* Array of XYZ values from colorimeter */
 static int create_ccmx(ccmx *p,
 char *desc,			/* General description (optional) */
 char *inst,			/* Instrument description to copy from */
-char *disp,			/* Display make and model (optional if tech) */
-char *tech,			/* Display technology description (optional if disp) */
+char *disp,			/* Display make and model (optional) */
+disptech dtech,		/* Display technology enum */
 int refrmode,		/* Display refresh mode, -1 = unknown, 0 = n, 1 = yes */
-int cbid,			/* Display type calibration base ID, 0 = unknown */
+int cbid,			/* Display type calibration base index, 0 = unknown */
 char *sel,			/* UI selector characters - NULL for none */
 char *refd,			/* Reference spectrometer description (optional) */
 int npat,			/* Number of samples in following arrays */
@@ -724,7 +698,7 @@ ccmx *new_ccmx(void) {
 		return NULL;
 
 	p->refrmode = -1;
-	p->cbid = 0;
+	p->cc_cbid = 0;
 
 	/* Init method pointers */
 	p->del            = del_ccmx;
diff --git a/xicc/ccmx.h b/xicc/ccmx.h
index b0ca01a..48598fc 100644
--- a/xicc/ccmx.h
+++ b/xicc/ccmx.h
@@ -33,13 +33,13 @@ struct _ccmx {
 	void (*del)(struct _ccmx *p);
 
 	/* Set the contents of the ccmx. return nz on error. */
-	int (*set_ccmx)(struct _ccmx *p, char *desc, char *inst, char *disp, char *tech,
+	int (*set_ccmx)(struct _ccmx *p, char *desc, char *inst, char *disp, disptech dtech,
 	                int refrmode, int cbid, char *sel, char *refd, double mtx[3][3]);	
 
 	/* Create a ccmx from measurements. return nz on error. */
-	int (*create_ccmx)(struct _ccmx *p, char *desc, char *inst, char *disp, char *tech,
+	int (*create_ccmx)(struct _ccmx *p, char *desc, char *inst, char *disp, disptech dtech,
 	               int refrmode, int cbid, char *sel, char *refd,
-	               int nsamples, double refs[][3], double cols[][3]);	
+	               int nsamples, double (*refs)[3], double (*cols)[3]);
 
 	/* write to a CGATS .ccmx file */
 	int (*write_ccmx)(struct _ccmx *p, char *filename);
@@ -64,10 +64,11 @@ struct _ccmx {
 	char *desc;		/* Desciption (optional) */
 	char *inst;		/* Name of colorimeter instrument */
 	char *disp;		/* Name of display (optional if tech) */
-	char *tech;		/* Technology (CRT, LCD + backlight type etc.) (optional if disp) */
-	int cbid;		/* Calibration display type base ID, 0 if not known */
+	disptech dtech;	/* Display Technology enumeration (optional if disp) */
+	char *tech;		/* Technology string (Looked up from dtech enum) */
+	int cc_cbid;	/* Calibration display type base ID required, 0 if not known */
 	int refrmode;	/* Refresh mode, -1 if unknown, 0 of no, 1 if yes */
-	char *sel;			/* Optional UI selector characters. May be NULL */
+	char *sel;		/* Optional UI selector characters. May be NULL */
 	char *ref;		/* Name of spectrometer instrument (optional) */
 	double matrix[3][3];	/* Transform matrix */
 	double av_err;			/* Average error of fit */
diff --git a/xicc/ccss.c b/xicc/ccss.c
index 6b5ebeb..663a789 100644
--- a/xicc/ccss.c
+++ b/xicc/ccss.c
@@ -35,6 +35,7 @@
 #endif
 #include "cgats.h"
 #include "xspect.h"
+#include "disptechs.h"
 #include "ccss.h"
 
 #ifdef NT       /* You'd think there might be some standards.... */
@@ -88,8 +89,9 @@ cgats **pocg		/* return CGATS structure */
 
 	if (p->disp)
 		ocg->add_kword(ocg, 0, "DISPLAY",p->disp, NULL);
-	if (p->tech)
-		ocg->add_kword(ocg, 0, "TECHNOLOGY", p->tech,NULL);
+
+	ocg->add_kword(ocg, 0, "TECHNOLOGY", disptech_get_id(p->dtech)->strid,NULL);
+
 	if (p->disp == NULL && p->tech == NULL) {
 		sprintf(p->err, "write_ccss: ccss doesn't contain display or techology strings");
 		ocg->del(ocg);
@@ -251,12 +253,10 @@ cgats *icg		/* input cgats structure */
 
 	if (icg->ntables == 0 || icg->t[0].tt != tt_other || icg->t[0].oi != 0) {
 		sprintf(p->err, "read_ccss: Input file isn't a CCSS format file");
-		icg->del(icg);
 		return 1;
 	}
 	if (icg->ntables != 1) {
 		sprintf(p->err, "Input file doesn't contain exactly one table");
-		icg->del(icg);
 		return 1;
 	}
 
@@ -265,21 +265,18 @@ cgats *icg		/* input cgats structure */
 	if ((ti = icg->find_kword(icg, 0, "DESCRIPTOR")) >= 0) {
 		if ((p->desc = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccss: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
 	if ((ti = icg->find_kword(icg, 0, "ORIGINATOR")) >= 0) {
 		if ((p->orig = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccss: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
 	if ((ti = icg->find_kword(icg, 0, "CREATED")) >= 0) {
 		if ((p->crdate = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccss: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
@@ -287,20 +284,19 @@ cgats *icg		/* input cgats structure */
 	if ((ti = icg->find_kword(icg, 0, "DISPLAY")) >= 0) {
 		if ((p->disp = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccss: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
 	if ((ti = icg->find_kword(icg, 0, "TECHNOLOGY")) >= 0) {
 		if ((p->tech = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccss: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
+		/* Get disptech enum from standard TECHNOLOGY string */
+		p->dtech = disptech_get_strid(p->tech)->dtech;
 	}
 	if (p->disp == NULL && p->tech == NULL) {
 		sprintf(p->err, "read_ccss: Input file doesn't contain keyword DISPLAY or TECHNOLOGY");
-		icg->del(icg);
 		return 1;
 	}
 	if ((ti = icg->find_kword(icg, 0, "DISPLAY_TYPE_REFRESH")) >= 0) {
@@ -313,7 +309,6 @@ cgats *icg		/* input cgats structure */
 	if ((ti = icg->find_kword(icg, 0, "UI_SELECTORS")) >= 0) {
 		if ((p->sel = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccss: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
@@ -321,26 +316,22 @@ cgats *icg		/* input cgats structure */
 	if ((ti = icg->find_kword(icg, 0, "REFERENCE")) >= 0) {
 		if ((p->ref = strdup(icg->t[0].kdata[ti])) == NULL) {
 			sprintf(p->err, "read_ccss: malloc failed");
-			icg->del(icg);
 			return 2;
 		}
 	}
 
 	if ((ii = icg->find_kword(icg, 0, "SPECTRAL_BANDS")) < 0) {
 		sprintf(p->err,"Input file doesn't contain keyword SPECTRAL_BANDS");
-		icg->del(icg);
 		return 1;
 	}
 	sp.spec_n = atoi(icg->t[0].kdata[ii]);
 	if ((ii = icg->find_kword(icg, 0, "SPECTRAL_START_NM")) < 0) {
 		sprintf(p->err,"Input file doesn't contain keyword SPECTRAL_START_NM");
-		icg->del(icg);
 		return 1;
 	}
 	sp.spec_wl_short = atof(icg->t[0].kdata[ii]);
 	if ((ii = icg->find_kword(icg, 0, "SPECTRAL_END_NM")) < 0) {
 		sprintf(p->err,"Input file doesn't contain keyword SPECTRAL_END_NM");
-		icg->del(icg);
 		return 1;
 	}
 	sp.spec_wl_long = atof(icg->t[0].kdata[ii]);
@@ -364,7 +355,6 @@ cgats *icg		/* input cgats structure */
 
 		if ((spi[j] = icg->find_field(icg, 0, buf)) < 0) {
 			sprintf(p->err,"Input file doesn't contain field %s",buf);
-			icg->del(icg);
 			return 1;
 		}
 	}
@@ -372,7 +362,6 @@ cgats *icg		/* input cgats structure */
 	if ((p->no_samp = icg->t[0].nsets) < 3) {
 		sprintf(p->err, "Input file doesn't contain at least three spectral samples");
 		p->no_samp = 0;
-		icg->del(icg);		/* Clean up */
 		return 1;
 	}
 
@@ -380,7 +369,6 @@ cgats *icg		/* input cgats structure */
 	if ((p->samples = (xspect *)malloc(sizeof(xspect) * p->no_samp)) == NULL) {
 		strcpy(p->err, "Malloc failed!");
 		p->no_samp = 0;
-		icg->del(icg);		/* Clean up */
 		return 2;
 	}
 
@@ -478,8 +466,8 @@ static int set_ccss(ccss *p,
 char *orig,			/* Originator (May be NULL) */
 char *crdate,		/* Creation date in ctime() format (May be NULL) */
 char *desc,			/* General description (optional) */
-char *disp,			/* Display make and model (optional if tech) */
-char *tech,			/* Display technology description (optional if disp) */
+char *disp,			/* Display make and model (optional) */
+disptech dtech,		/* Display technology enum */
 int refrmode,		/* Display refresh mode, -1 = unknown, 0 = n, 1 = yes */
 char *sel,			/* UI selector characters - NULL for none */
 char *refd,			/* Reference spectrometer description (optional) */
@@ -513,12 +501,7 @@ int no_samp			/* Number of spectral samples */
 			return 2;
 		}
 	}
-	if (tech != NULL) {
-		if ((p->tech = strdup(tech)) == NULL) {
-			sprintf(p->err, "set_ccss: malloc tech failed");
-			return 2;
-		}
-	}
+	p->dtech = dtech;
 	p->refrmode = refrmode;
 	if (sel != NULL) {
 		if ((p->sel = strdup(sel)) == NULL) {
diff --git a/xicc/ccss.h b/xicc/ccss.h
index 0d2e128..924f045 100644
--- a/xicc/ccss.h
+++ b/xicc/ccss.h
@@ -33,7 +33,7 @@ struct _ccss {
 	/* Set the contents of the ccss. return nz on error. */
 	/* (Makes copies of all parameters) */
 	int (*set_ccss)(struct _ccss *p, char *orig, char *cdate,
-	                char *desc, char *disp, char *tech, int refrmode, char *sel,
+	                char *desc, char *disp, disptech dtech, int refrmode, char *sel,
 	                char *ref, xspect *samples, int no_samp);	
 
 	/* write to a CGATS .ccss file */
@@ -58,7 +58,8 @@ struct _ccss {
 	char *crdate;		/* Creation date (in ctime() format). May be NULL */
 	char *desc;			/* General Description (optional) */
 	char *disp;			/* Description of the display (Manfrr and Model No) (optional if tech) */
-	char *tech;			/* Technology (CRT, LCD + backlight type etc.) (optional if disp) */
+	disptech dtech;		/* Display Technology enumeration (optional if disp) */
+	char *tech;			/* Technology string (Looked up from dtech enum) */
 	int refrmode;		/* Refresh mode, -1 if unknown, 0 of no, 1 if yes */
 	char *sel;			/* Optional UI selector characters. May be NULL */
 	char *ref;			/* Name of reference spectrometer instrument (optional) */
diff --git a/xicc/ccttest.c b/xicc/ccttest.c
index 9662678..28a3dc9 100644
--- a/xicc/ccttest.c
+++ b/xicc/ccttest.c
@@ -23,6 +23,7 @@
 #include "xspect.h"
 #include "numlib.h"
 #include "plot.h"
+#include "ui.h"
 
 #define PLANKIAN
 #define XRES 500
diff --git a/xicc/cgatsplot.c b/xicc/cgatsplot.c
index 9fadf89..ccaa28a 100644
--- a/xicc/cgatsplot.c
+++ b/xicc/cgatsplot.c
@@ -23,9 +23,10 @@
 #include "numlib.h"
 #include "icc.h"
 #include "cgats.h"
-#include "plot.h"
 #include "xcolorants.h"
 #include "sort.h"
+#include "plot.h"
+#include "ui.h"
 
 void usage(void) {
 	fprintf(stderr,"Simple 2D plot of CGATS .ti3 data\n");
diff --git a/xicc/cv.c b/xicc/cv.c
index 5271bc4..7cd7927 100644
--- a/xicc/cv.c
+++ b/xicc/cv.c
@@ -17,6 +17,7 @@
 #include <math.h>
 #include <numlib.h>
 #include "plot.h"
+#include "ui.h"
 
 void usage(void);
 
diff --git a/xicc/cvtest.c b/xicc/cvtest.c
index 6d24572..a9fa59f 100644
--- a/xicc/cvtest.c
+++ b/xicc/cvtest.c
@@ -27,6 +27,7 @@
 #endif
 #include "numlib.h"
 #include "plot.h"
+#include "ui.h"
 
 #define MAX_PARM 40		/* Make > SHAPE_ORDS */
 
diff --git a/xicc/extractttag.c b/xicc/extractttag.c
index 79b4085..2c26bda 100644
--- a/xicc/extractttag.c
+++ b/xicc/extractttag.c
@@ -1,6 +1,7 @@
 
 /* 
  * Extract a CGATS file from an ICC profile tag.
+ * (Can also extract a tag of unknown format as a binary lump).
  *
  * Author:  Graeme W. Gill
  * Date:    2008/5/18
@@ -14,6 +15,9 @@
 
 /*
  * TTBD:
+ *
+ * Should uncompress ZXML type tag using zlib.
+ *
  */
 
 
@@ -28,6 +32,7 @@
 #include "numlib.h"
 #include "icc.h"
 #include "xicc.h"
+#include "ui.h"
 
 #define MXTGNMS 30
 
@@ -60,6 +65,7 @@ main(int argc, char *argv[]) {
 	icc *icco;
 	icTagSignature sig;
 	icmText *ro;
+	icmUnknown *uro;
 	icmFile *ifp, *ofp;
 	int verb = 0;
 	int  size = 0;
@@ -135,11 +141,14 @@ main(int argc, char *argv[]) {
 
 	sig = str2tag(tag_name);
 
-	if ((ro = (icmText *)icco->read_tag(icco, sig)) == NULL)
+	if ((ro = (icmText *)icco->read_tag_any(icco, sig)) == NULL) {
 		error("%d, %s",icco->errc, icco->err);
+	}
 
-	if (ro->ttype != icSigTextType) {
-		error("Tag isn't TextType");
+	if (ro->ttype == icmSigUnknownType) {
+		uro = (icmUnknown *)ro;
+	} else if (ro->ttype != icSigTextType) {
+		error("Tag isn't TextType or UnknownType");
 	}
 
 	if (docal) {
@@ -187,8 +196,14 @@ main(int argc, char *argv[]) {
 			error("unable to open output file '%s'",out_name);
 		}
 	
-		if (ofp->write(ofp, ro->data, 1, ro->size-1) != (ro->size-1)) {
-			error("writing to file '%s' failed",out_name);
+		if (ro->ttype == icmSigUnknownType) {
+			if (ofp->write(ofp, uro->data, 1, uro->size) != (uro->size)) {
+				error("writing to file '%s' failed",out_name);
+			}
+		} else {
+			if (ofp->write(ofp, ro->data, 1, ro->size-1) != (ro->size-1)) {
+				error("writing to file '%s' failed",out_name);
+			}
 		}
 	
 		if (ofp->del(ofp) != 0) {
diff --git a/xicc/fakeCMY.c b/xicc/fakeCMY.c
index 4387cd8..b27f2d1 100644
--- a/xicc/fakeCMY.c
+++ b/xicc/fakeCMY.c
@@ -31,6 +31,7 @@
 #include "xicc.h"
 #include "cgats.h"
 #include "targen.h"
+#include "ui.h"
 
 #define EXTRA_NEUTRAL 50		/* Crude way of increasing weighting */
 
diff --git a/xicc/fbview.c b/xicc/fbview.c
index 53da99e..07a9170 100644
--- a/xicc/fbview.c
+++ b/xicc/fbview.c
@@ -25,6 +25,7 @@
 #include "aconfig.h"
 #include "numlib.h"
 #include "icc.h"
+#include "vrml.h"
 
 #define TRES 43
 
@@ -124,7 +125,7 @@ main(
 	strcpy(out_name, in_name);
 	if ((xl = strrchr(out_name, '.')) == NULL)	/* Figure where extention is */
 		xl = out_name + strlen(out_name);
-	strcpy(xl,".wrl");
+	xl[0] = '\000';			/* Remove extension */
 
 	/* Open up the file for reading */
 	if ((rd_fp = new_icmFileStd_name(in_name,"r")) == NULL)
@@ -145,18 +146,8 @@ main(
 		double aerr = 0.0;
 		double nsamps = 0.0;
 		icmLuBase *luo1, *luo2;
-		FILE *wrl;
-		struct {
-			double x, y, z;
-			double wx, wy, wz;
-			double r, g, b;
-		} axes[5] = {
-			{ 0, 0,  50-GAMUT_LCENT,  2, 2, 100,  .7, .7, .7 },	/* L axis */
-			{ 50, 0,  0-GAMUT_LCENT,  100, 2, 2,   1,  0,  0 },	/* +a (red) axis */
-			{ 0, -50, 0-GAMUT_LCENT,  2, 100, 2,   0,  0,  1 },	/* -b (blue) axis */
-			{ -50, 0, 0-GAMUT_LCENT,  100, 2, 2,   0,  1,  0 },	/* -a (green) axis */
-			{ 0,  50, 0-GAMUT_LCENT,  2, 100, 2,   1,  1,  0 },	/* +b (yellow) axis */
-		};
+		int doaxes = 1;
+		vrml *wrl;
 		int i, j;
 	
 
@@ -180,53 +171,11 @@ main(
 				error ("%d, %s",rd_icco->errc, rd_icco->err);
 		}
 
-		if ((wrl = fopen(out_name,"w")) == NULL) {
-			fprintf(stderr,"Error opening output file '%s'\n",out_name);
-			return 2;
-		}
-
-		/* Spit out a VRML 2 Object surface of gamut */
-
-		fprintf(wrl,"#VRML V2.0 utf8\n");
-		fprintf(wrl,"\n");
-		fprintf(wrl,"# Created by the Argyll CMS\n");
-		fprintf(wrl,"Transform {\n");
-		fprintf(wrl,"children [\n");
-		fprintf(wrl,"	NavigationInfo {\n");
-		fprintf(wrl,"		type \"EXAMINE\"        # It's an object we examine\n");
-		fprintf(wrl,"	} # We'll add our own light\n");
-		fprintf(wrl,"\n");
-		fprintf(wrl,"    DirectionalLight {\n");
-		fprintf(wrl,"        direction 0 0 -1      # Light illuminating the scene\n");
-		fprintf(wrl,"        direction 0 -1 0      # Light illuminating the scene\n");
-		fprintf(wrl,"    }\n");
-		fprintf(wrl,"\n");
-		fprintf(wrl,"    Viewpoint {\n");
-		fprintf(wrl,"        position 0 0 340      # Position we view from\n");
-		fprintf(wrl,"    }\n");
-		fprintf(wrl,"\n");
-		if (doaxes != 0) {
-			fprintf(wrl,"# Lab axes as boxes:\n");
-			for (i = 0; i < 5; i++) {
-				fprintf(wrl,"Transform { translation %f %f %f\n", axes[i].x, axes[i].y, axes[i].z);
-				fprintf(wrl,"\tchildren [\n");
-				fprintf(wrl,"\t\tShape{\n");
-				fprintf(wrl,"\t\t\tgeometry Box { size %f %f %f }\n",
-				                  axes[i].wx, axes[i].wy, axes[i].wz);
-				fprintf(wrl,"\t\t\tappearance Appearance { material Material ");
-				fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[i].r, axes[i].g, axes[i].b);
-				fprintf(wrl,"\t\t}\n");
-				fprintf(wrl,"\t]\n");
-				fprintf(wrl,"}\n");
-			}
-			fprintf(wrl,"\n");
+		if ((wrl = new_vrml(out_name, doaxes, vrml_lab)) == NULL) {
+			error("new_vrml for '%s%s' failed\n",out_name,vrml_ext());
 		}
 
-		fprintf(wrl,"\n");
-		fprintf(wrl,"Shape {\n");
-		fprintf(wrl,"  geometry IndexedLineSet { \n");
-		fprintf(wrl,"    coord Coordinate { \n");
-		fprintf(wrl,"	   point [\n");
+		wrl->start_line_set(wrl, 0);
 
 		i = 0;
 //		for (co[0] = 0; co[0] < TRES; co[0]++) {
@@ -269,33 +218,21 @@ main(
 					aerr += absd;
 
 					if (absd > 3.0) {
-						fprintf(wrl,"%f %f %f,\n",in[1], in[2], in[0]-GAMUT_LCENT);
-						fprintf(wrl,"%f %f %f,\n",check[1], check[2], check[0]-GAMUT_LCENT);
+						wrl->add_vertex(wrl, 0, in);
+						wrl->add_vertex(wrl, 0, check);
 						i++;
 					}
 				}
 			}
 //		}
 
-		fprintf(wrl,"      ]\n");
-		fprintf(wrl,"    }\n");
-		fprintf(wrl,"  coordIndex [\n");
+		wrl->make_lines(wrl, 0, 2);
 
-		for (j = 0; j < i; j++) {
-			fprintf(wrl,"%d, %d, -1,\n", j * 2, j * 2 + 1);
-		}
-		fprintf(wrl,"    ]\n");
-		fprintf(wrl,"  }\n");
-		fprintf(wrl,"} # end shape\n");
-
-		fprintf(wrl,"\n");
-		fprintf(wrl,"  ] # end of children for world\n");
-		fprintf(wrl,"}\n");
-	
-		if (fclose(wrl) != 0) {
-			fprintf(stderr,"Error closing output file '%s'\n",out_name);
+		if (wrl->flush(wrl) != 0) {
+			fprintf(stderr,"Error writint output file '%s%s'\n",out_name,vrml_ext());
 			return 2;
 		}
+		wrl->del(wrl);
 
 		printf("bwd to fwd check complete, peak err = %f, avg err = %f\n",merr,aerr/nsamps);
 
diff --git a/xicc/iccgamut.c b/xicc/iccgamut.c
index bb97f26..a77653f 100644
--- a/xicc/iccgamut.c
+++ b/xicc/iccgamut.c
@@ -41,6 +41,8 @@
 #include "xicc.h"
 #include "gamut.h"
 #include "counters.h"
+#include "vrml.h"
+#include "ui.h"
 
 static void diag_gamut(icxLuBase *p, double detail, int doaxes,
                        double tlimit, double klimit, char *outname);
@@ -54,9 +56,10 @@ void usage(char *diag) {
 		fprintf(stderr,"Diagnostic: %s\n",diag);
 	fprintf(stderr," -v            Verbose\n");
 	fprintf(stderr," -d sres       Surface resolution details 1.0 - 50.0\n");
-	fprintf(stderr," -w            emit VRML .wrl file as well as CGATS .gam file\n");
-	fprintf(stderr," -n            Don't add VRML axes or white/black point\n");
-	fprintf(stderr," -k            Add VRML markers for prim. & sec. \"cusp\" points\n");
+	fprintf(stderr," -w            emit %s %s file as well as CGATS .gam file\n",vrml_format(),vrml_ext());
+	fprintf(stderr," -n            Don't add %s axes or white/black point\n",vrml_format());
+	fprintf(stderr," -k            Add %s markers for prim. & sec. \"cusp\" points\n",vrml_format());
+	fprintf(stderr,"               (Set env. ARGYLL_3D_DISP_FORMAT to VRML, X3D or X3DOM to change format)\n");
 	fprintf(stderr," -f function   f = forward*, b = backwards\n");
 	fprintf(stderr," -i intent     p = perceptual, r = relative colorimetric,\n");
 	fprintf(stderr,"               s = saturation, a = absolute (default), d = profile default\n");
@@ -479,20 +482,20 @@ main(int argc, char *argv[]) {
 	if (special) {
 		if (func != icmFwd)
 			error("Must be forward direction for special plot");
-		strcpy(xl,".wrl");
+		xl[0] = '\000';			/* remove extension */
 		diag_gamut(luo, gamres, doaxes, tlimit/100.0, klimit/100.0, out_name); 
 	} else {
 		/* Creat a gamut surface */
 		if ((gam = luo->get_gamut(luo, gamres)) == NULL)
 			error ("%d, %s",xicco->errc, xicco->err);
 
-		if (gam->write_gam(gam,out_name))
+		if (gam->write_gam(gam, out_name))
 			error ("write gamut failed on '%s'",out_name);
 
 		if (vrml) {
-			strcpy(xl,".wrl");
+			xl[0] = '\000';			/* remove extension */
 			if (gam->write_vrml(gam,out_name, doaxes, docusps))
-				error ("write vrml failed on '%s'",out_name);
+				error ("write vrml failed on '%s%s'",out_name, vrml_ext());
 		}
 
 		if (verb) {
@@ -524,21 +527,10 @@ double detail,		/* Gamut resolution detail */
 int doaxes,			/* Do Lab axes */
 double tlimit,		/* Total ink limit */
 double klimit,		/* K ink limit */
-char *outname		/* Output VRML file */
+char *outname		/* Output VRML/X3D file (no extension) */
 ) {
 	int i, j;
-	FILE *wrl;
-	struct {
-		double x, y, z;
-		double wx, wy, wz;
-		double r, g, b;
-	} axes[5] = {
-		{ 0, 0,  50-GAMUT_LCENT,  2, 2, 100,  .7, .7, .7 },	/* L axis */
-		{ 50, 0,  0-GAMUT_LCENT,  100, 2, 2,   1,  0,  0 },	/* +a (red) axis */
-		{ 0, -50, 0-GAMUT_LCENT,  2, 100, 2,   0,  0,  1 },	/* -b (blue) axis */
-		{ -50, 0, 0-GAMUT_LCENT,  100, 2, 2,   0,  1,  0 },	/* -a (green) axis */
-		{ 0,  50, 0-GAMUT_LCENT,  2, 100, 2,   1,  1,  0 },	/* +b (yellow) axis */
-	};
+	vrml *wrl;
 	int vix;						/* Vertex index */
 	DCOUNT(coa, MXDI, p->inputChan, 0, 0, 2);
 
@@ -589,55 +581,10 @@ char *outname		/* Output VRML file */
 	if (res < 2)
 		res = 2;
 
-	if ((wrl = fopen(outname,"w")) == NULL)
-		error("Error opening wrl output file '%s'",outname);
-
-	/* Spit out a VRML 2 Object surface of gamut */
-	fprintf(wrl,"#VRML V2.0 utf8\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"# Created by the Argyll CMS\n");
-	fprintf(wrl,"Transform {\n");
-  	fprintf(wrl,"children [\n");
-    fprintf(wrl,"	NavigationInfo {\n");
-	fprintf(wrl,"		type \"EXAMINE\"        # It's an object we examine\n");
-	fprintf(wrl,"	} # We'll add our own light\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"    DirectionalLight {\n");
-	fprintf(wrl,"        direction 0 0 -1      # Light illuminating the scene\n");
-	fprintf(wrl,"        direction 0 -1 0      # Light illuminating the scene\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"    Viewpoint {\n");
-	fprintf(wrl,"        position 0 0 340      # Position we view from\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"\n");
-	if (doaxes != 0) {
-		fprintf(wrl,"# Lab axes as boxes:\n");
-		for (i = 0; i < 5; i++) {
-			fprintf(wrl,"Transform { translation %f %f %f\n", axes[i].x, axes[i].y, axes[i].z);
-			fprintf(wrl,"\tchildren [\n");
-			fprintf(wrl,"\t\tShape{\n");
-			fprintf(wrl,"\t\t\tgeometry Box { size %f %f %f }\n",
-			                  axes[i].wx, axes[i].wy, axes[i].wz);
-			fprintf(wrl,"\t\t\tappearance Appearance { material Material ");
-			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[i].r, axes[i].g, axes[i].b);
-			fprintf(wrl,"\t\t}\n");
-			fprintf(wrl,"\t]\n");
-			fprintf(wrl,"}\n");
-		}
-		fprintf(wrl,"\n");
-	}
-	fprintf(wrl,"    Transform {\n");
-	fprintf(wrl,"      translation 0 0 0\n");
-	fprintf(wrl,"      children [\n");
-	fprintf(wrl,"		Shape { \n");
-	fprintf(wrl,"		    geometry IndexedFaceSet {\n");
-	fprintf(wrl,"				solid FALSE\n");		/* Don't back face cull */
-	fprintf(wrl,"				convex TRUE\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"		        coord Coordinate { \n");
-	fprintf(wrl,"		            point [			# Verticy coordinates\n");
+	if ((wrl = new_vrml(outname, doaxes, vrml_lab)) == NULL)
+		error("Error creating wrl output '%s%s'",outname,vrml_ext());
 
+	wrl->start_line_set(wrl, 0);		/* Start set 0 */
 
 	/* Itterate over all the faces in the device space */
 	/* generating the vertx positions. */
@@ -645,7 +592,7 @@ char *outname		/* Output VRML file */
 	vix = 0;
 	while(!DC_DONE(coa)) {
 		int e, m1, m2;
-		double in[MXDI];
+		double in[MXDI], xb, yb;
 		double inl[MXDI];
 		double out[3];
 		double sum;
@@ -666,111 +613,28 @@ char *outname		/* Output VRML file */
 
 				/* Scan over 2D device space face */
 				for (x = 0; x < res; x++) {				/* step over surface */
-					in[m1] = x/(res - 1.0);
+					xb = in[m1] = x/(res - 1.0);
 					for (y = 0; y < res; y++) {
-						in[m2] = y/(res - 1.0);
+						double rgb[3], rgb2[3];
+						int v0, v1, v2, v3;
+						int ix[4];
 
-						for (sum = 0.0, e = 0; e < p->inputChan; e++) {
+						yb = in[m2] = y/(res - 1.0);
+
+						/* Check ink limit */
+						for (sum = 0.0, e = 0; e < p->inputChan; e++)
 							sum += inl[e] = in[e];
-						}
 						if (sum >= tlimit) {
 							for (e = 0; e < p->inputChan; e++)
 								inl[e] *= tlimit/sum;
 						}
 						if (p->inputChan >= 3 && inl[3] >= klimit)
 							inl[3] = klimit;
-						p->lookup(p, out, inl);
-						fprintf(wrl,"%f %f %f,\n",out[1], out[2], out[0]-50.0);
-						vix++;
-					}
-				}
-			}
-		}
-		/* Increment index within block */
-		DC_INC(coa);
-	}
-
-	fprintf(wrl,"					]\n");
-	fprintf(wrl,"		        }\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"		        coordIndex [ 		# Indexes of poligon Verticies \n");
-
-	/* Itterate over all the faces in the device space */
-	/* generating the quadrilateral indexes. */
-	DC_INIT(coa);
-	vix = 0;
-	while(!DC_DONE(coa)) {
-		int e, m1, m2;
-		double in[MXDI];
-
-		/* Scan only device surface */
-		for (m1 = 0; m1 < p->inputChan; m1++) {
-			if (coa[m1] != 0)
-				continue;
-
-			for (m2 = m1 + 1; m2 < p->inputChan; m2++) {
-				int x, y;
-
-				if (coa[m2] != 0)
-					continue;
-
-				for (e = 0; e < p->inputChan; e++)
-					in[e] = (double)coa[e];		/* Base value */
-
-				/* Scan over 2D device space face */
-				/* Only output quads under the total ink limit */
-				/* Scan over 2D device space face */
-				for (x = 0; x < res; x++) {				/* step over surface */
-					for (y = 0; y < res; y++) {
-						if (x < (res-1) && y < (res-1)) {
-							fprintf(wrl,"%d, %d, %d, %d, -1\n", 
-							vix, vix + 1, vix + 1 + res, vix + res);
-						}
-						vix++;
-					}
-				}
-			}
-		}
-		/* Increment index within block */
-		DC_INC(coa);
-	}
-
-	fprintf(wrl,"				]\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"				colorPerVertex TRUE\n");
-	fprintf(wrl,"		        color Color {\n");
-	fprintf(wrl,"		            color [			# RGB colors of each vertex\n");
-
-	/* Itterate over all the faces in the device space */
-	/* generating the vertx colors. */
-	DC_INIT(coa);
-	vix = 0;
-	while(!DC_DONE(coa)) {
-		int e, m1, m2;
-		double in[MXDI];
-
-		/* Scan only device surface */
-		for (m1 = 0; m1 < p->inputChan; m1++) {
-			if (coa[m1] != 0)
-				continue;
-
-			for (m2 = m1 + 1; m2 < p->inputChan; m2++) {
-				int x, y;
-
-				if (coa[m2] != 0)
-					continue;
-
-				for (e = 0; e < p->inputChan; e++)
-					in[e] = (double)coa[e];		/* Base value */
 
-				/* Scan over 2D device space face */
-				for (x = 0; x < res; x++) {				/* step over surface */
-					double xb = x/(res - 1.0);
-					for (y = 0; y < res; y++) {
-						int v0, v1, v2, v3;
-						double yb = y/(res - 1.0);
-						double rgb[3];
+						/* Lookup L*a*b* value */
+						p->lookup(p, out, inl);
 
+						/* Compute color */
 						for (v0 = 0, e = 0; e < p->inputChan; e++)
 							v0 |= coa[e] ? (1 << e) : 0;		/* Binary index */
 
@@ -785,7 +649,20 @@ char *outname		/* Output VRML file */
 							       + (1.0 - yb) *        xb  * col[v3][j]
 							       +        yb  *        xb  * col[v2][j];
 						}
-						fprintf(wrl,"%f %f %f,\n",rgb[1], rgb[2], rgb[0]);
+						/* re-order the color */
+						rgb2[0] = rgb[1];
+						rgb2[1] = rgb[2];
+						rgb2[2] = rgb[0];
+						wrl->add_col_vertex(wrl, 0, out, rgb2);
+
+						/* Add the quad vertexes */
+						if (x < (res-1) && y < (res-1)) {
+							ix[0] = vix;
+							ix[1] = vix + 1;
+							ix[2] = vix + 1 + res;
+							ix[3] = vix + res;
+							wrl->add_quad(wrl, 0, ix);
+						}
 						vix++;
 					}
 				}
@@ -795,25 +672,8 @@ char *outname		/* Output VRML file */
 		DC_INC(coa);
 	}
 
-	fprintf(wrl,"					] \n");
-	fprintf(wrl,"		        }\n");
-	fprintf(wrl,"		    }\n");
-	fprintf(wrl,"		    appearance Appearance { \n");
-	fprintf(wrl,"		        material Material {\n");
-	fprintf(wrl,"					transparency 0.0\n");
-	fprintf(wrl,"					ambientIntensity 0.3\n");
-	fprintf(wrl,"					shininess 0.5\n");
-	fprintf(wrl,"				}\n");
-	fprintf(wrl,"		    }\n");
-	fprintf(wrl,"		}	# end Shape\n");
-	fprintf(wrl,"      ]\n");
-	fprintf(wrl,"    }\n");
-
-	fprintf(wrl,"\n");
-	fprintf(wrl,"  ] # end of children for world\n");
-	fprintf(wrl,"}\n");
-
-	if (fclose(wrl) != 0)
-		error("Error closing output file '%s'",outname);
+	wrl->make_quads_vc(wrl, 0, 0.0);	/* Make set 0 with color per vertex and 0 transparence */
+
+	wrl->del(wrl);
 }
 
diff --git a/xicc/icheck.c b/xicc/icheck.c
index ed3c6e5..5f61d9a 100644
--- a/xicc/icheck.c
+++ b/xicc/icheck.c
@@ -12,6 +12,12 @@
  * Please refer to License.txt file for details.
  */
 
+/*
+	Estimate PCS->device interpolation error by comparing
+	the interpolated PCS value in the cell to the PCS computed
+	from the corners of the cell.
+*/
+
 /* TTBD:
  *
  */
@@ -26,23 +32,19 @@
 #include "copyright.h"
 #include "aconfig.h"
 #include "icc.h"
+#include "vrml.h"
 
 void usage(void) {
 	fprintf(stderr,"Check PCS->Device Interpolation faults of ICC file, Version %s\n",ARGYLL_VERSION_STR);
 	fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
 	fprintf(stderr,"usage: icheck [-v] [-w] infile\n");
 	fprintf(stderr," -v        verbose\n");
-	fprintf(stderr," -w        create VRML visualisation\n");
-	fprintf(stderr," -x        Use VRML axies\n");
+	fprintf(stderr," -a        Show all interpolation errors, not just worst\n");
+	fprintf(stderr," -w        create %s visualisation\n",vrml_format());
+	fprintf(stderr," -x        Use %s axies\n",vrml_format());
 	exit(1);
 }
 
-FILE *start_vrml(char *name, int doaxes);
-void start_line_set(FILE *wrl);
-void add_vertex(FILE *wrl, double pp[3]);
-void make_lines(FILE *wrl, int ppset);
-void end_vrml(FILE *wrl);
-
 int
 main(
 	int argc,
@@ -50,6 +52,7 @@ main(
 ) {
 	int fa,nfa;				/* argument we're looking at */
 	int verb = 0;
+	int doall = 0;
 	int dovrml = 0;
 	int doaxes = 0;
 	char in_name[100];
@@ -65,7 +68,7 @@ main(
 	icColorSpaceSignature ins, outs;	/* Type of input and output spaces */
 	int inn;							/* Number of input chanels */
 	icmLuAlgType alg;
-	FILE *wrl = NULL;
+	vrml *wrl = NULL;
 	
 	error_program = argv[0];
 
@@ -93,7 +96,7 @@ main(
 			if (argv[fa][1] == 'v' || argv[fa][1] == 'V') {
 				verb = 1;
 			}
-			/* VRML */
+			/* VRML/X3D */
 			else if (argv[fa][1] == 'w' || argv[fa][1] == 'W') {
 				dovrml = 1;
 			}
@@ -101,6 +104,9 @@ main(
 			else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') {
 				doaxes = 1;
 			}
+			else if (argv[fa][1] == 'a') {
+				doall = 1;
+			}
 			else if (argv[fa][1] == '?')
 				usage();
 			else 
@@ -116,7 +122,7 @@ main(
 	strcpy(out_name, in_name);
 	if ((xl = strrchr(out_name, '.')) == NULL)	/* Figure where extention is */
 		xl = out_name + strlen(out_name);
-	strcpy(xl,".wrl");
+	xl[0] = '\000';				/* Remove extension */
 
 	/* Open up the file for reading */
 	if ((rd_fp = new_icmFileStd_name(in_name,"r")) == NULL)
@@ -158,8 +164,10 @@ main(
 	gres = lluob->lut->clutPoints;
 
 	if (dovrml) {
-		wrl = start_vrml(out_name, doaxes);
-		start_line_set(wrl);
+		wrl = new_vrml(out_name, doaxes, vrml_lab);
+		if (wrl == NULL)
+			error("new_vrml failed for '%s%s'",out_name,vrml_ext());
+		wrl->start_line_set(wrl, 0);
 	}
 
 	{
@@ -188,9 +196,9 @@ main(
 					/* For each corner of the PCS grid based at the current point, */
 					/* average the PCS and Device values */
 					m = 0;
-					for (cc[2] = 0; cc[2] < 2; cc[2]++, m++) {
+					for (cc[2] = 0; cc[2] < 2; cc[2]++) {
 						for (cc[1] = 0; cc[1] < 2; cc[1]++) {
-							for (cc[0] = 0; cc[0] < 2; cc[0]++) {
+							for (cc[0] = 0; cc[0] < 2; cc[0]++, m++) {
 								double dev[MAX_CHAN];
 
 								pcs[m][0] = (co[0] + cc[0])/(gres - 1.0);
@@ -231,7 +239,7 @@ main(
 					for (k = 0; k < inn; k++)
 						adev[k] /= 8.0;
 	
-					/* Compute worst case distance of PCS corners to average PCS */
+					/* Compute worst case distance of PCS corners to PCS of average dev */
 					wpcsd = 0.0;
 					for (m = 0; m < 8; m++) {
 						double ss;
@@ -243,7 +251,10 @@ main(
 						if (ss > wpcsd)
 							wpcsd = ss;
 					}
-					wpcsd *= 0.75;		/* Set threshold at 75% of most distant corner */
+					wpcsd *= 0.15;		/* Set threshold at 75% of most distant corner */
+
+//printf("~1 wpcsd = %f\n",wpcsd);
+
 					/* Set a worst case */
 					if (wpcsd < 1.0)
 						wpcsd = 1.0;
@@ -277,14 +288,14 @@ main(
 					if (ier > merr)
 						merr = ier;
 
-					if (ier > wpcsd) {
+					if (doall || ier > wpcsd) {
 						tcount++;
 
 						printf("ier = %f, Dev = %f %f %f %f\n",
 						       ier, adev[0], adev[1], adev[2], adev[3]);
 						if (dovrml) {
-							add_vertex(wrl, apcs);
-							add_vertex(wrl, check);
+							wrl->add_vertex(wrl, 0, apcs);
+							wrl->add_vertex(wrl, 0, check);
 						}
 					}
 
@@ -299,8 +310,8 @@ main(
 		}
 
 		if (dovrml) {
-			make_lines(wrl, 2);
-			end_vrml(wrl);
+			wrl->make_lines(wrl, 0, 2);
+			wrl->del(wrl);
 		}
 	
 		aerr /= ccount;
@@ -319,214 +330,3 @@ main(
 	return 0;
 }
 
-/* ------------------------------------------------ */
-/* Some simple functions to do basix VRML work */
-
-#define GAMUT_LCENT 50.0
-static int npoints = 0;
-static int paloc = 0;
-static struct { double pp[3]; } *pary;
-
-static void Lab2RGB(double *out, double *in);
-
-FILE *start_vrml(char *name, int doaxes) {
-	FILE *wrl;
-	struct {
-		double x, y, z;
-		double wx, wy, wz;
-		double r, g, b;
-	} axes[5] = {
-		{ 0, 0,  50-GAMUT_LCENT,  2, 2, 100,  .7, .7, .7 },	/* L axis */
-		{ 50, 0,  0-GAMUT_LCENT,  100, 2, 2,   1,  0,  0 },	/* +a (red) axis */
-		{ 0, -50, 0-GAMUT_LCENT,  2, 100, 2,   0,  0,  1 },	/* -b (blue) axis */
-		{ -50, 0, 0-GAMUT_LCENT,  100, 2, 2,   0,  1,  0 },	/* -a (green) axis */
-		{ 0,  50, 0-GAMUT_LCENT,  2, 100, 2,   1,  1,  0 },	/* +b (yellow) axis */
-	};
-	int i;
-	
-	if ((wrl = fopen(name,"w")) == NULL)
-		error("Error opening VRML file '%s'\n",name);
-
-	npoints = 0;
-
-	fprintf(wrl,"#VRML V2.0 utf8\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"# Created by the Argyll CMS\n");
-	fprintf(wrl,"Transform {\n");
-	fprintf(wrl,"children [\n");
-	fprintf(wrl,"	NavigationInfo {\n");
-	fprintf(wrl,"		type \"EXAMINE\"        # It's an object we examine\n");
-	fprintf(wrl,"	} # We'll add our own light\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"    DirectionalLight {\n");
-	fprintf(wrl,"        direction 0 0 -1      # Light illuminating the scene\n");
-	fprintf(wrl,"        direction 0 -1 0      # Light illuminating the scene\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"    Viewpoint {\n");
-	fprintf(wrl,"        position 0 0 340      # Position we view from\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"\n");
-	if (doaxes != 0) {
-		fprintf(wrl,"# Lab axes as boxes:\n");
-		for (i = 0; i < 5; i++) {
-			fprintf(wrl,"Transform { translation %f %f %f\n", axes[i].x, axes[i].y, axes[i].z);
-			fprintf(wrl,"\tchildren [\n");
-			fprintf(wrl,"\t\tShape{\n");
-			fprintf(wrl,"\t\t\tgeometry Box { size %f %f %f }\n",
-			                  axes[i].wx, axes[i].wy, axes[i].wz);
-			fprintf(wrl,"\t\t\tappearance Appearance { material Material ");
-			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[i].r, axes[i].g, axes[i].b);
-			fprintf(wrl,"\t\t}\n");
-			fprintf(wrl,"\t]\n");
-			fprintf(wrl,"}\n");
-		}
-		fprintf(wrl,"\n");
-	}
-
-	return wrl;
-}
-
-void
-start_line_set(FILE *wrl) {
-
-	fprintf(wrl,"\n");
-	fprintf(wrl,"Shape {\n");
-	fprintf(wrl,"  geometry IndexedLineSet { \n");
-	fprintf(wrl,"    coord Coordinate { \n");
-	fprintf(wrl,"	   point [\n");
-}
-
-void add_vertex(FILE *wrl, double pp[3]) {
-
-	fprintf(wrl,"%f %f %f,\n",pp[1], pp[2], pp[0]-GAMUT_LCENT);
-	
-	if (paloc < (npoints+1)) {
-		paloc = (paloc + 10) * 2;
-		if (pary == NULL)
-			pary = malloc(paloc * 3 * sizeof(double));
-		else
-			pary = realloc(pary, paloc * 3 * sizeof(double));
-
-		if (pary == NULL)
-			error ("Malloc failed");
-	}
-	pary[npoints].pp[0] = pp[0];
-	pary[npoints].pp[1] = pp[1];
-	pary[npoints].pp[2] = pp[2];
-	npoints++;
-}
-
-
-void make_lines(FILE *wrl, int ppset) {
-	int i, j;
-
-	fprintf(wrl,"      ]\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"  coordIndex [\n");
-
-	for (i = 0; i < npoints;) {
-		for (j = 0; j < ppset; j++, i++) {
-			fprintf(wrl,"%d, ", i);
-		}
-		fprintf(wrl,"-1,\n");
-	}
-	fprintf(wrl,"    ]\n");
-
-	/* Color */
-	fprintf(wrl,"            colorPerVertex TRUE\n");
-	fprintf(wrl,"            color Color {\n");
-	fprintf(wrl,"              color [			# RGB colors of each vertex\n");
-
-	for (i = 0; i < npoints; i++) {
-		double rgb[3], Lab[3];
-		Lab[0] = pary[i].pp[0];
-		Lab[1] = pary[i].pp[1];
-		Lab[2] = pary[i].pp[2];
-		Lab2RGB(rgb, Lab);
-		fprintf(wrl,"                %f %f %f,\n", rgb[0], rgb[1], rgb[2]);
-	}
-	fprintf(wrl,"              ] \n");
-	fprintf(wrl,"            }\n");
-	/* End color */
-
-	fprintf(wrl,"  }\n");
-	fprintf(wrl,"} # end shape\n");
-
-}
-
-void end_vrml(FILE *wrl) {
-
-	fprintf(wrl,"\n");
-	fprintf(wrl,"  ] # end of children for world\n");
-	fprintf(wrl,"}\n");
-
-	if (fclose(wrl) != 0)
-		error("Error closing VRML file\n");
-}
-
-
-/* Convert a gamut Lab value to an RGB value for display purposes */
-static void
-Lab2RGB(double *out, double *in) {
-	double L = in[0], a = in[1], b = in[2];
-	double x,y,z,fx,fy,fz;
-	double R, G, B;
-
-	/* Scale so that black is visible */
-	L = L * (100 - 40.0)/100.0 + 40.0;
-
-	/* First convert to XYZ using D50 white point */
-	if (L > 8.0) {
-		fy = (L + 16.0)/116.0;
-		y = pow(fy,3.0);
-	} else {
-		y = L/903.2963058;
-		fy = 7.787036979 * y + 16.0/116.0;
-	}
-
-	fx = a/500.0 + fy;
-	if (fx > 24.0/116.0)
-		x = pow(fx,3.0);
-	else
-		x = (fx - 16.0/116.0)/7.787036979;
-
-	fz = fy - b/200.0;
-	if (fz > 24.0/116.0)
-		z = pow(fz,3.0);
-	else
-		z = (fz - 16.0/116.0)/7.787036979;
-
-	x *= 0.9642;	/* Multiply by white point, D50 */
-	y *= 1.0;
-	z *= 0.8249;
-
-	/* Now convert to sRGB values */
-	R = x * 3.2410  + y * -1.5374 + z * -0.4986;
-	G = x * -0.9692 + y * 1.8760  + z * 0.0416;
-	B = x * 0.0556  + y * -0.2040 + z * 1.0570;
-
-	if (R < 0.0)
-		R = 0.0;
-	else if (R > 1.0)
-		R = 1.0;
-
-	if (G < 0.0)
-		G = 0.0;
-	else if (G > 1.0)
-		G = 1.0;
-
-	if (B < 0.0)
-		B = 0.0;
-	else if (B > 1.0)
-		B = 1.0;
-
-	R = pow(R, 1.0/2.2);
-	G = pow(G, 1.0/2.2);
-	B = pow(B, 1.0/2.2);
-
-	out[0] = R;
-	out[1] = G;
-	out[2] = B;
-}
-
diff --git a/xicc/monctest.c b/xicc/monctest.c
index 8630099..a97b1b5 100644
--- a/xicc/monctest.c
+++ b/xicc/monctest.c
@@ -26,8 +26,9 @@
 #include "copyright.h"
 #include "aconfig.h"
 #include "numlib.h"
-#include "plot.h"
 #include "moncurve.h"
+#include "plot.h"
+#include "ui.h"
 
 double lin(double x, double xa[], double ya[], int n);
 void usage(void);
@@ -141,17 +142,24 @@ int main() {
 
 			/* Create X values */
 			for (i = 0; i < pnts; i++)
-				xa[i] = i/(pnts-1.0);
+				xa[i] = pow(i/(pnts-1.0), 1.0);
+
+			for (i = 0; i < pnts; i++)
+				ya[i] = pow(xa[i], ex);
+
+		} else if (n == 1) {				/* inverse exponential function aproximation */
+			double ex = 1.0/2.4;
+			pnts = MAX_PNTS;
+
+			printf("Trial %d, no points = %d, exponential %f\n",n,pnts,ex);
+
+			/* Create X values */
+			for (i = 0; i < pnts; i++)
+				xa[i] = pow(i/(pnts-1.0), 3.0);
 
 			for (i = 0; i < pnts; i++)
 				ya[i] = pow(xa[i], ex);
 
-		} else if (n < (TSETS+1))	/* Standard versions */ {
-			pnts = t1p[n-1];
-			for (i = 0; i < pnts; i++) {
-				xa[i] = t1xa[n-1][i];
-				ya[i] = t1ya[n-1][i];
-			}
 #endif
 
 		} else {	/* Random versions */
diff --git a/xicc/mpp.c b/xicc/mpp.c
index 02b3019..86b6c30 100644
--- a/xicc/mpp.c
+++ b/xicc/mpp.c
@@ -90,9 +90,17 @@
 						/* accuracy worse!) */
 
 						/* Transfer curve parameter (wiggle) minimisation weight */
-#define TRANS_BASE  0.2	/* 0 & 1 harmonic parameter weight */
-#define TRANS_HBASE 0.8	/* 2nd harmonic and above base parameter weight */
-#define SHAPE_PMW 0.2	/* Shape parameter (wiggle) minimisation weight */
+//#define TRANS_BASE  0.2	/* 0 & 1 harmonic parameter weight */
+//#define TRANS_HBASE 0.8	/* 2nd harmonic and above base parameter weight */
+#define TRANS_HW01		  0.01		/* 0 & 1 harmonic weights */
+#define TRANS_HBREAK	  3			/* Harmonic that has HWBR */
+//#define TRANS_HWBR        0.5		/* Base weight of harmonics HBREAK up */
+//#define TRANS_HWINC       0.5		/* Increase in weight for each harmonic above HWBR */
+#define TRANS_HWBR        0.5		/* Base weight of harmonics HBREAK up */
+#define TRANS_HWINC       0.5		/* Increase in weight for each harmonic above HWBR */
+
+//#define SHAPE_PMW 0.2	/* Shape parameter (wiggle) minimisation weight */
+#define SHAPE_PMW 10.0	/* Shape parameter (wiggle) minimisation weight */
 #define COMB_PMW 0.008	/* Primary combination anchor point distance weight */
 
 #define verbo stdout
@@ -1916,7 +1924,20 @@ static double efunc2(void *adata, double pv[]) {
 		double w;
 		for (k = 0; k < p->cord; k++) {
 			i = m * p->cord + k;
+
+#ifdef NEVER
 			w = (k < 2) ? TRANS_BASE : k * TRANS_HBASE; 	/* Increase weight with harmonics */
+#else
+			/* Weigh to suppress ripples */
+			if (k <= 1) {						/* Use TRANS_HW01 */
+				w = TRANS_HW01;
+			} else if (k <= TRANS_HBREAK) {	/* Blend from TRANS_HW01 to TRANS_HWBR */
+				double bl = (k - 1.0)/(TRANS_HBREAK - 1.0);
+				w = (1.0 - bl) * TRANS_HW01 + bl * TRANS_HWBR;
+			} else {				/* Use TRANS_HWBR */
+				w = TRANS_HWBR + (k-TRANS_HBREAK) * TRANS_HWINC;
+			}
+#endif
 			smv += w * pv[i] * pv[i];
 		}
 	}
@@ -2224,7 +2245,19 @@ for (m = 0; m < p->n; m++) {
 		double w;
 		for (k = 0; k < p->cord; k++) {
 			i = m * p->cord + k;
+#ifdef NEVER
 			w = (k < 2) ? TRANS_BASE : k * TRANS_HBASE; 	/* Increase weight with harmonics */
+#else
+			/* Weigh to suppress ripples */
+			if (k <= 1) {						/* Use TRANS_HW01 */
+				w = TRANS_HW01;
+			} else if (k <= TRANS_HBREAK) {	/* Blend from TRANS_HW01 to TRANS_HWBR */
+				double bl = (k - 1.0)/(TRANS_HBREAK - 1.0);
+				w = (1.0 - bl) * TRANS_HW01 + bl * TRANS_HWBR;
+			} else {				/* Use TRANS_HWBR */
+				w = TRANS_HWBR + (k-TRANS_HBREAK) * TRANS_HWINC;
+			}
+#endif
 			dv[i] += w * tt * pv[i];				/* Del in rv due to del in pv */
 			smv += w * pv[i] * pv[i];
 		}
@@ -3597,24 +3630,24 @@ static int create(
 	switch (quality) {
 		case 0:				/* Low */
 			useshape = 0;
-			mxtcord = 3;
-			maxit = 2;
+			mxtcord = 8;
+			maxit = 3;
 			break;
 		case 1:
 		default:			/* Medium */
 			useshape = 1;
-			mxtcord = 4;
-			maxit = 2;
+			mxtcord = 10;
+			maxit = 4;
 			break;
 		case 2:				/* High */
 			useshape = 1;
-			mxtcord = 5;
-			maxit = 3;
+			mxtcord = 14;
+			maxit = 5;
 			break;
 		case 3:				/* Ultra high */
 			useshape = 1;
-			mxtcord = 10;	/* Is more actually better ? */
-			maxit = 8;
+			mxtcord = 20;	/* Is more actually better ? */
+			maxit = 10;
 			break;
 		case 99:				/* Special, simple model */
 			useshape = 0;
diff --git a/xicc/mpp.h b/xicc/mpp.h
index 826c4de..26fa3a8 100644
--- a/xicc/mpp.h
+++ b/xicc/mpp.h
@@ -27,7 +27,7 @@
 /* ------------------------------------------------------------------------------ */
 
 #define MPP_MXINKS 8		/* Would like to be ICX_MXINKS but need more dynamic allocation */
-#define MPP_MXTCORD 10		/* Maximum shaper harmonic orders to use */
+#define MPP_MXTCORD 20		/* Maximum shaper harmonic orders to use */
 #define MPP_MXCCOMB (1 << MPP_MXINKS)	/* Maximum number of primary combinations */
 #define MPP_MXPARMS (MPP_MXINKS * MPP_MXTCORD + (MPP_MXINKS * MPP_MXCCOMB/2) + MPP_MXCCOMB)
 							/* Maximum total parameters for a band */
diff --git a/xicc/mpplu.c b/xicc/mpplu.c
index a2bb9ad..11ddb76 100644
--- a/xicc/mpplu.c
+++ b/xicc/mpplu.c
@@ -26,6 +26,8 @@
 #include "numlib.h"
 #include "xicc.h"
 #include "counters.h"
+#include "vrml.h"
+#include "ui.h"
 
 void usage(void) {
 	fprintf(stderr,"Translate colors through an MPP profile, V1.00\n");
@@ -41,13 +43,13 @@ void usage(void) {
 	fprintf(stderr,"            1931_2 (def), 1964_10, S&B 1955_2, shaw, J&V 1978_2\n");
 	fprintf(stderr," -u         Use Fluorescent Whitening Agent compensation\n");
 	fprintf(stderr," -g         Create gamut output\n");
-	fprintf(stderr," -w         Create gamut VRML as well\n");
-	fprintf(stderr," -n         Don't add VRML axes\n");
+	fprintf(stderr," -w         Create gamut %s as well\n",vrml_format());
+	fprintf(stderr," -n         Don't add %s axes\n",vrml_format());
 	fprintf(stderr," -a n       Gamut transparency level\n");
 	fprintf(stderr," -d n       Gamut surface detail level\n");
 	fprintf(stderr," -t num     Invoke debugging test code \"num\" 1..n\n");
 	fprintf(stderr,"            1 - check partial derivative for device input\n");
-	fprintf(stderr,"            2 - create overlap diagnostic VRML gamut surface\n");
+	fprintf(stderr,"            2 - create overlap diagnostic %s gamut surface\n",vrml_format());
 	fprintf(stderr,"\n");
 	fprintf(stderr,"    The colors to be translated should be fed into stdin,\n");
 	fprintf(stderr,"    one input color per line, white space separated.\n");
@@ -67,8 +69,8 @@ main(int argc, char *argv[]) {
 	int verb = 0;
 	int test = 0;			/* special test code */
 	int dogam = 0;			/* Create gamut */
-	int dowrl = 0;			/* Create VRML gamut */
-	int doaxes = 1;			/* Create VRML axes */
+	int dowrl = 0;			/* Create VRML/X3D gamut */
+	int doaxes = 1;			/* Create VRML/X3D axes */
 	double trans = 0.0;		/* Transparency */
 	double gamres = 0.0;	/* Gamut resolution */
 	int repYxy = 0;			/* Report Yxy */
@@ -250,13 +252,13 @@ main(int argc, char *argv[]) {
 			else if (argv[fa][1] == 'g' || argv[fa][1] == 'G')
 				dogam = 1;
 
-			/* VRML plot */
+			/* VRML/X3D plot */
 			else if (argv[fa][1] == 'w' || argv[fa][1] == 'W') {
 				dogam = 1;
 				dowrl = 1;
 			}
 
-			/* No VRML axes */
+			/* No VRML/X3D axes */
 			else if (argv[fa][1] == 'n' || argv[fa][1] == 'N') {
 				doaxes = 0;
 			}
@@ -415,8 +417,7 @@ main(int argc, char *argv[]) {
 			strcpy(gam_name, prof_name);
 			if ((xl = strrchr(gam_name, '.')) == NULL)	/* Figure where extention is */
 				xl = gam_name + strlen(gam_name);
-	
-			strcpy(xl,".wrl");
+			xl[0] = '\000';				/* Remove extension */
 			diag_gamut(mppo, gamres, doaxes, trans, gam_name);
 
 		} else {
@@ -434,15 +435,14 @@ main(int argc, char *argv[]) {
 		strcpy(gam_name, prof_name);
 		if ((xl = strrchr(gam_name, '.')) == NULL)	/* Figure where extention is */
 			xl = gam_name + strlen(gam_name);
-	
 		strcpy(xl,".gam");
-		if (gam->write_gam(gam,gam_name))
+		if (gam->write_gam(gam, gam_name))
 			error ("write gamut failed on '%s'",gam_name);
 	
 		if (dowrl) {
-			strcpy(xl,".wrl");
+			xl[0] = '\000';
 			if (gam->write_vrml(gam,gam_name, doaxes, docusps))
-				error ("write vrml failed on '%s'",gam_name);
+				error ("write vrml failed on '%s%s'",gam_name,vrml_ext());
 		}
 
 		gam->del(gam);
@@ -602,21 +602,10 @@ mpp *p,				/* This */
 double detail,		/* Gamut resolution detail */
 int doaxes,
 double trans,		/* Transparency */
-char *outname		/* Output VRML file */
+char *outname		/* Output VRML/X3D file (no extension) */
 ) {
 	int i, j;
-	FILE *wrl;
-	struct {
-		double x, y, z;
-		double wx, wy, wz;
-		double r, g, b;
-	} axes[5] = {
-		{ 0, 0,  50-GAMUT_LCENT,  2, 2, 100,  .7, .7, .7 },	/* L axis */
-		{ 50, 0,  0-GAMUT_LCENT,  100, 2, 2,   1,  0,  0 },	/* +a (red) axis */
-		{ 0, -50, 0-GAMUT_LCENT,  2, 100, 2,   0,  0,  1 },	/* -b (blue) axis */
-		{ -50, 0, 0-GAMUT_LCENT,  100, 2, 2,   0,  1,  0 },	/* -a (green) axis */
-		{ 0,  50, 0-GAMUT_LCENT,  2, 100, 2,   1,  1,  0 },	/* +b (yellow) axis */
-	};
+	vrml *wrl;
 	int vix;						/* Vertex index */
 	DCOUNT(coa, MAX_CHAN, p->n, 0, 0, 2);
 	double col[MPP_MXCCOMB][3];		/* Color asigned to each major vertex */
@@ -661,55 +650,10 @@ char *outname		/* Output VRML file */
 	if (res < 2)
 		res = 2;
 
-	if ((wrl = fopen(outname,"w")) == NULL)
-		error("Error opening wrl output file '%s'",outname);
-
-	/* Spit out a VRML 2 Object surface of gamut */
-	fprintf(wrl,"#VRML V2.0 utf8\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"# Created by the Argyll CMS\n");
-	fprintf(wrl,"Transform {\n");
-  	fprintf(wrl,"children [\n");
-    fprintf(wrl,"	NavigationInfo {\n");
-	fprintf(wrl,"		type \"EXAMINE\"        # It's an object we examine\n");
-	fprintf(wrl,"	} # We'll add our own light\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"    DirectionalLight {\n");
-	fprintf(wrl,"        direction 0 0 -1      # Light illuminating the scene\n");
-	fprintf(wrl,"        direction 0 -1 0      # Light illuminating the scene\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"    Viewpoint {\n");
-	fprintf(wrl,"        position 0 0 340      # Position we view from\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"\n");
-	if (doaxes != 0) {
-		fprintf(wrl,"# Lab axes as boxes:\n");
-		for (i = 0; i < 5; i++) {
-			fprintf(wrl,"Transform { translation %f %f %f\n", axes[i].x, axes[i].y, axes[i].z);
-			fprintf(wrl,"\tchildren [\n");
-			fprintf(wrl,"\t\tShape{\n");
-			fprintf(wrl,"\t\t\tgeometry Box { size %f %f %f }\n",
-			                  axes[i].wx, axes[i].wy, axes[i].wz);
-			fprintf(wrl,"\t\t\tappearance Appearance { material Material ");
-			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[i].r, axes[i].g, axes[i].b);
-			fprintf(wrl,"\t\t}\n");
-			fprintf(wrl,"\t]\n");
-			fprintf(wrl,"}\n");
-		}
-		fprintf(wrl,"\n");
-	}
-	fprintf(wrl,"    Transform {\n");
-	fprintf(wrl,"      translation 0 0 0\n");
-	fprintf(wrl,"      children [\n");
-	fprintf(wrl,"		Shape { \n");
-	fprintf(wrl,"		    geometry IndexedFaceSet {\n");
-	fprintf(wrl,"				solid FALSE\n");		/* Don't back face cull */
-	fprintf(wrl,"				convex TRUE\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"		        coord Coordinate { \n");
-	fprintf(wrl,"		            point [			# Verticy coordinates\n");
+	if ((wrl = new_vrml(outname, doaxes, vrml_lab)) == NULL)
+		error("new_vrml faile for file '%s%s'",outname,vrml_ext());
 
+	wrl->start_line_set(wrl, 0);
 
 	/* Itterate over all the faces in the device space */
 	/* generating the vertx positions. */
@@ -720,6 +664,7 @@ char *outname		/* Output VRML file */
 		double in[MAX_CHAN];
 		double out[3];
 		double sum;
+		double xb, yb;
 
 		/* Scan only device surface */
 		for (m1 = 0; m1 < p->n; m1++) {
@@ -737,57 +682,32 @@ char *outname		/* Output VRML file */
 
 				/* Scan over 2D device space face */
 				for (x = 0; x < res; x++) {				/* step over surface */
-					in[m1] = x/(res - 1.0);
+					xb = in[m1] = x/(res - 1.0);
 					for (y = 0; y < res; y++) {
-						in[m2] = y/(res - 1.0);
+						int v0, v1, v2, v3;
+						double rgb[3];
+						yb = in[m2] = y/(res - 1.0);
 
+						/* Lookup PCS value */
 						p->lookup(p, out, in);
-						fprintf(wrl,"%f %f %f,\n",out[1], out[2], out[0]-50.0);
-						vix++;
-					}
-				}
-			}
-		}
-		/* Increment index within block */
-		DC_INC(coa);
-	}
-
-	fprintf(wrl,"					]\n");
-	fprintf(wrl,"		        }\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"		        coordIndex [ 		# Indexes of poligon Verticies \n");
 
-	/* Itterate over all the faces in the device space */
-	/* generating the quadrilateral indexes. */
-	DC_INIT(coa);
-	vix = 0;
-	while(!DC_DONE(coa)) {
-		int e, m1, m2;
-		double in[MAX_CHAN];
-		double sum;
-
-		/* Scan only device surface */
-		for (m1 = 0; m1 < p->n; m1++) {
-			if (coa[m1] != 0)
-				continue;
-
-			for (m2 = m1 + 1; m2 < p->n; m2++) {
-				int x, y;
-
-				if (coa[m2] != 0)
-					continue;
-
-				for (sum = 0.0, e = 0; e < p->n; e++)
-					in[e] = (double)coa[e];		/* Base value */
+						/* Create a color */
+						for (v0 = 0, e = 0; e < p->n; e++)
+							v0 |= coa[e] ? (1 << e) : 0;		/* Binary index */
 
-				/* Scan over 2D device space face */
-				for (x = 0; x < res; x++) {				/* step over surface */
-					for (y = 0; y < res; y++) {
+						v1 = v0 | (1 << m2);				/* Y offset */
+						v2 = v0 | (1 << m2) | (1 << m1);	/* X+Y offset */
+						v3 = v0 | (1 << m1);				/* Y offset */
 
-						if (x < (res-1) && y < (res-1)) {
-							fprintf(wrl,"%d, %d, %d, %d, -1\n", 
-							vix, vix + 1, vix + 1 + res, vix + res);
+						/* Linear interp between the main verticies */
+						for (j = 0; j < 3; j++) {
+							rgb[j] = (1.0 - yb) * (1.0 - xb) * col[v0][j]
+							       +        yb  * (1.0 - xb) * col[v1][j]
+							       + (1.0 - yb) *        xb  * col[v3][j]
+							       +        yb  *        xb  * col[v2][j];
 						}
+
+						wrl->add_col_vertex(wrl, 0, out, rgb);
 						vix++;
 					}
 				}
@@ -797,14 +717,8 @@ char *outname		/* Output VRML file */
 		DC_INC(coa);
 	}
 
-	fprintf(wrl,"				]\n");
-	fprintf(wrl,"\n");
-	fprintf(wrl,"				colorPerVertex TRUE\n");
-	fprintf(wrl,"		        color Color {\n");
-	fprintf(wrl,"		            color [			# RGB colors of each vertex\n");
-
 	/* Itterate over all the faces in the device space */
-	/* generating the vertx colors. */
+	/* generating the quadrilateral indexes. */
 	DC_INIT(coa);
 	vix = 0;
 	while(!DC_DONE(coa)) {
@@ -828,27 +742,16 @@ char *outname		/* Output VRML file */
 
 				/* Scan over 2D device space face */
 				for (x = 0; x < res; x++) {				/* step over surface */
-					double xb = x/(res - 1.0);
 					for (y = 0; y < res; y++) {
-						int v0, v1, v2, v3;
-						double yb = y/(res - 1.0);
-						double rgb[3];
-
-						for (v0 = 0, e = 0; e < p->n; e++)
-							v0 |= coa[e] ? (1 << e) : 0;		/* Binary index */
-
-						v1 = v0 | (1 << m2);				/* Y offset */
-						v2 = v0 | (1 << m2) | (1 << m1);	/* X+Y offset */
-						v3 = v0 | (1 << m1);				/* Y offset */
 
-						/* Linear interp between the main verticies */
-						for (j = 0; j < 3; j++) {
-							rgb[j] = (1.0 - yb) * (1.0 - xb) * col[v0][j]
-							       +        yb  * (1.0 - xb) * col[v1][j]
-							       + (1.0 - yb) *        xb  * col[v3][j]
-							       +        yb  *        xb  * col[v2][j];
+						if (x < (res-1) && y < (res-1)) {
+							int ix[4];
+							ix[0] = vix;
+							ix[1] = vix + 1;
+							ix[2] = vix + 1 + res;
+							ix[3] = vix + res;
+							wrl->add_quad(wrl, 0, ix);
 						}
-						fprintf(wrl,"%f %f %f,\n",rgb[1], rgb[2], rgb[0]);
 						vix++;
 					}
 				}
@@ -858,26 +761,12 @@ char *outname		/* Output VRML file */
 		DC_INC(coa);
 	}
 
-	fprintf(wrl,"					] \n");
-	fprintf(wrl,"		        }\n");
-	fprintf(wrl,"		    }\n");
-	fprintf(wrl,"		    appearance Appearance { \n");
-	fprintf(wrl,"		        material Material {\n");
-	fprintf(wrl,"					transparency %f\n",trans);
-	fprintf(wrl,"					ambientIntensity 0.3\n");
-	fprintf(wrl,"					shininess 0.5\n");
-	fprintf(wrl,"				}\n");
-	fprintf(wrl,"		    }\n");
-	fprintf(wrl,"		}	# end Shape\n");
-	fprintf(wrl,"      ]\n");
-	fprintf(wrl,"    }\n");
-
-	fprintf(wrl,"\n");
-	fprintf(wrl,"  ] # end of children for world\n");
-	fprintf(wrl,"}\n");
-
-	if (fclose(wrl) != 0)
-		error("Error closing output file '%s'",outname);
+	wrl->make_quads_vc(wrl, 0, trans);
+
+	if (wrl->flush(wrl) != 0)
+		error("Error closing output file '%s%s'",outname,vrml_ext());
+
+	wrl->del(wrl);
 }
 
 /* -------------------------------------------- */
diff --git a/xicc/revfix.c b/xicc/revfix.c
index ee9662a..2add264 100644
--- a/xicc/revfix.c
+++ b/xicc/revfix.c
@@ -34,6 +34,7 @@
 #include "aconfig.h"
 #include "numlib.h"
 #include "xicc.h"
+#include "ui.h"
 
 #define USE_CAM_CLIP_OPT		/* Clip in CAM Jab space rather than Lab */
 #undef DEBUG		/* Print each value changed */
diff --git a/xicc/specplot.c b/xicc/specplot.c
index 8e52726..9de4645 100644
--- a/xicc/specplot.c
+++ b/xicc/specplot.c
@@ -24,6 +24,7 @@
 #include "xspect.h"
 #include "numlib.h"
 #include "plot.h"
+#include "ui.h"
 
 #define PLANKIAN
 #define XRES 500
diff --git a/xicc/specsubsamp.c b/xicc/specsubsamp.c
index b13ad94..4163ee2 100644
--- a/xicc/specsubsamp.c
+++ b/xicc/specsubsamp.c
@@ -21,7 +21,7 @@
 #include <math.h>
 #include "xspect.h"
 #include "numlib.h"
-
+#include "ui.h"
 
 void usage(void) {
 	fprintf(stderr,"Downsample spectral data\n");
diff --git a/xicc/spectest.c b/xicc/spectest.c
index 1c92df4..c340e98 100644
--- a/xicc/spectest.c
+++ b/xicc/spectest.c
@@ -34,6 +34,7 @@
 #ifdef DOPLOT
 #include "plot.h"
 #endif
+#include "ui.h"
 
 
 /* Spectrolino filter "D65" illuminant */
diff --git a/xicc/spectest2.c b/xicc/spectest2.c
index e9c1754..3fd574f 100644
--- a/xicc/spectest2.c
+++ b/xicc/spectest2.c
@@ -32,6 +32,7 @@
 #ifdef DOPLOT
 #include "plot.h"
 #endif
+#include "ui.h"
 
 
 /* Normal 'A' spectra, then UV filtered version */
diff --git a/xicc/tiffgamut.c b/xicc/tiffgamut.c
index 9423eae..9ddb62b 100644
--- a/xicc/tiffgamut.c
+++ b/xicc/tiffgamut.c
@@ -22,6 +22,9 @@
  *      Need to cope with profile not having black point.
  *
  *		How to cope with no profile, therefore no white or black point ?
+ *
+ *		Should we have a median filter option, to ignore small groups
+ *		of extreme pixel values, rathar than total small numbers using -f ?
  */
 
 
@@ -42,6 +45,8 @@
 #include "gamut.h"
 #include "xicc.h"
 #include "sort.h"
+#include "vrml.h"
+#include "ui.h"
 
 #undef NOCAMGAM_CLIP		/* No clip to CAM gamut before CAM lookup */
 #undef DEBUG				/* Dump filter cell contents */
@@ -58,14 +63,15 @@ void del_filter();
 
 void usage(void) {
 	int i;
-	fprintf(stderr,"Create VRML image of the gamut surface of a TIFF or JPEG, Version %s\n",ARGYLL_VERSION_STR);
+	fprintf(stderr,"Create gamut surface of a TIFF or JPEG, Version %s\n",ARGYLL_VERSION_STR);
 	fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
 	fprintf(stderr,"usage: tiffgamut [-v level] [profile.icm | embedded.tif/jpg] infile1.tif/jpg [infile2.tif/jpg ...] \n");
 	fprintf(stderr," -v            Verbose\n");
 	fprintf(stderr," -d sres       Surface resolution details 1.0 - 50.0\n");
-	fprintf(stderr," -w            emit VRML .wrl file as well as CGATS .gam file\n");
-	fprintf(stderr," -n            Don't add VRML axes or white/black point\n");
-	fprintf(stderr," -k            Add markers for prim. & sec. \"cusp\" points\n");
+	fprintf(stderr," -w            emit %s %s file as well as CGATS .gam file\n",vrml_format(),vrml_ext());
+	fprintf(stderr," -n            Don't add %s axes or white/black point\n",vrml_format());
+	fprintf(stderr," -k            Add %s markers for prim. & sec. \"cusp\" points\n",vrml_format());
+	fprintf(stderr,"               (set env. ARGYLL_3D_DISP_FORMAT to VRML, X3D or X3DOM to change format)\n");
 	fprintf(stderr," -f perc       Filter by popularity, perc = percent to use\n");
 	fprintf(stderr," -i intent     p = perceptual, r = relative colorimetric,\n");
 	fprintf(stderr,"               s = saturation, a = absolute (default), d = profile default\n");
@@ -325,7 +331,7 @@ main(int argc, char *argv[]) {
 	int ffa, lfa;				/* First, last input file argument */
 	char prof_name[MAXNAMEL+1] = { '\000' };	/* ICC profile name, "" if none */
 	char in_name[MAXNAMEL+1];			/* TIFF input file */
-	char *xl = NULL, out_name[MAXNAMEL+4+1] = { '\000' };	/* VRML output file */
+	char *xl = NULL, out_name[MAXNAMEL+4+1] = { '\000' };	/* VRML/X3D output file */
 	int verb = 0;
 	int vrml = 0;
 	int doaxes = 1;
@@ -366,12 +372,12 @@ main(int argc, char *argv[]) {
 	uint16 samplesperpixel, bitspersample;
 	uint16 pconfig, photometric, pmtc;
 	tdata_t *inbuf;
-	int inbpix;                					/* Number of pixels in jpeg in buf */
+	int inbpix = 0;               					/* Number of pixels in jpeg in buf */
 	void (*cvt)(double *out, double *in) = NULL;	/* TIFF conversion function, NULL if none */
-	icColorSpaceSignature tcs;					/* TIFF colorspace */
-	uint16 extrasamples;						/* Extra "alpha" samples */
-	uint16 *extrainfo;							/* Info about extra samples */
-	int sign_mask;								/* Handling of encoding sign */
+	icColorSpaceSignature tcs = 0;				/* TIFF colorspace */
+	uint16 extrasamples = 0;					/* Extra "alpha" samples */
+	uint16 *extrainfo = NULL;					/* Info about extra samples */
+	int sign_mask = 0;							/* Handling of encoding sign */
 
 	/* JPEG */
 	jpegerrorinfo jpeg_rerr;
@@ -547,7 +553,7 @@ main(int argc, char *argv[]) {
 					usage();
 			}
 
-			/* VRML output */
+			/* VRML/X3D output */
 			else if (argv[fa][1] == 'w' || argv[fa][1] == 'W') {
 				vrml = 1;
 			}
@@ -1010,19 +1016,26 @@ main(int argc, char *argv[]) {
 				int i;
 				double in[MAX_CHAN], out[MAX_CHAN];
 				
+//printf("~1 location %d,%d\n",x,y);
 				if (bitspersample == 8) {
 					for (i = 0; i < samplesperpixel; i++) {
 						int v = ((unsigned char *)inbuf)[x * samplesperpixel + i];
+//printf("~1 v[%d] = %d\n",i,v);
 						if (sign_mask & (1 << i))		/* Treat input as signed */
 							v = (v & 0x80) ? v - 0x80 : v + 0x80;
+//printf("~1 signed v[%d] = %d\n",i,v);
 						in[i] = v/255.0;
+//printf("~1 in[%d] = %f\n",i,in[i]);
 					}
 				} else {
 					for (i = 0; i < samplesperpixel; i++) {
 						int v = ((unsigned short *)inbuf)[x * samplesperpixel + i];
+//printf("~1 v[%d] = %d\n",i,v);
 						if (sign_mask & (1 << i))		/* Treat input as signed */
 							v = (v & 0x8000) ? v - 0x8000 : v + 0x8000;
+//printf("~1 signed v[%d] = %d\n",i,v);
 						in[i] = v/65535.0;
+//printf("~1 in[%d] = %f\n",i,in[i]);
 					}
 				}
 				if (cvt != NULL) {	/* Undo TIFF encoding */
@@ -1054,10 +1067,14 @@ main(int argc, char *argv[]) {
 				}
 
 				for (i = 0; i < 3; i++) {
-					if (out[i] < apcsmin[i])
+					if (out[i] < apcsmin[i]) {
+//printf("~1 new min %f\n",out[i]);
 						apcsmin[i] = out[i];
-					if (out[i] > apcsmax[i])
+					}
+					if (out[i] > apcsmax[i]) {
+//printf("~1 new max %f\n",out[i]);
 						apcsmax[i] = out[i];
+					}
 				}
 				if (filter)
 					add_fpixel(out);
@@ -1117,14 +1134,13 @@ main(int argc, char *argv[]) {
 	if (verb)
 		printf("Output Gamut file '%s'\n",out_name);
 
-	/* Create the VRML file */
+	/* Create the VRML/X3D file */
 	if (gam->write_gam(gam,out_name))
 		error ("write gamut failed on '%s'",out_name);
 
 	if (vrml) {
-
-		strcpy(xl,".wrl");
-		printf("Output vrml file '%s'\n",out_name);
+		xl[0] = '\000';
+		printf("Output %s file '%s%s'\n",vrml_format(),out_name,vrml_ext());
 		if (gam->write_vrml(gam,out_name, doaxes, docusps))
 			error ("write vrml failed on '%s'",out_name);
 	}
diff --git a/xicc/tiffgmts.c b/xicc/tiffgmts.c
index 03f6363..c31f852 100644
--- a/xicc/tiffgmts.c
+++ b/xicc/tiffgmts.c
@@ -1,6 +1,6 @@
 
 /* 
- * Create a gamut mapping test set from a TIFF file.
+ * Create a gamut mapping locus set from a TIFF file.
  *
  * Author:  Graeme W. Gill
  * Date:    08/10/14
@@ -42,6 +42,7 @@
 #include "xicc.h"
 #include "vrml.h"
 #include "sort.h"
+#include "ui.h"
 
 #define DE_SPACE 3		/* Delta E of spacing for output points */
 #undef DEBUG_PLOT
@@ -52,8 +53,8 @@ void usage(void) {
 	fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
 	fprintf(stderr,"usage: tiffgmts [-v level] [profile.icm | embedded.tif] infile.tif\n");
 	fprintf(stderr," -v            Verbose\n");
-	fprintf(stderr," -w            emit VRML .wrl file as well as CGATS .ts file\n");
-	fprintf(stderr," -n            Don't add VRML axes or white/black point\n");
+	fprintf(stderr," -w            emit %s %s file as well as CGATS .ts file\n",vrml_format(),vrml_ext());
+	fprintf(stderr," -n            Don't add %s axes or white/black point\n",vrml_format());
 	fprintf(stderr," -i intent     p = perceptual, r = relative colorimetric,\n");
 	fprintf(stderr,"               s = saturation, a = absolute (default), d = profile default\n");
 //  fprintf(stderr,"               P = absolute perceptual, S = absolute saturation\n");
@@ -484,7 +485,7 @@ main(int argc, char *argv[]) {
 					usage();
 			}
 
-			/* VRML output */
+			/* VRML/X3D output */
 			else if (argv[fa][1] == 'w' || argv[fa][1] == 'W') {
 				dovrml = 1;
 			}
@@ -1003,14 +1004,14 @@ printf("~1 itter %d, alen = %f, minl = %f, maxl = %f\n",j,alen,minl,maxl);
 			    error("Write error : %s",pp->err);
 		}
 
-		/* Create the VRML file */
+		/* Create the VRML/X3D file */
 		if (dovrml) {
 			vrml *vv;
 			
-			strcpy(xl,".wrl");
-			printf("Output vrml file '%s'\n",out_name);
+			xl[0] = '\000';			/* remove extension */
+			printf("Output %s file '%s%s'\n",vrml_format(),out_name,vrml_ext());
 			if ((vv = new_vrml(out_name, doaxes, 0)) == NULL)
-				error ("Creating VRML object failed");
+				error ("Creating %s object %s%s failed",vrml_format(),out_name,vrml_ext());
 
 #ifdef NEVER
 			vv->start_line_set(vv);
@@ -1024,7 +1025,7 @@ printf("~1 itter %d, alen = %f, minl = %f, maxl = %f\n",j,alen,minl,maxl);
 			}
 #endif
 
-			vv->del(vv);
+			vv->del(vv);		/* Write file */
 		}
 		free(outp);
 	}
diff --git a/xicc/transplot.c b/xicc/transplot.c
index d4600d9..6473e2b 100644
--- a/xicc/transplot.c
+++ b/xicc/transplot.c
@@ -1,6 +1,5 @@
 
 /* 
- * International Color Consortium Format Library (icclib)
  * Check various aspects of RGB or CMYK device link, 
  * and RGB/CMYK profile transfer characteristics.
  *
@@ -29,6 +28,7 @@
 #include "icc.h"
 #include "numlib.h"
 #include "plot.h"
+#include "ui.h"
 
 void usage(void) {
 	fprintf(stderr,"Check CMYK/RGB/PCS->PCS/CMYK/RGB transfer response\n");
@@ -321,8 +321,12 @@ main(
 		}
 		if (labout)
 			do_plot6(xx,y0,y1,NULL,NULL,y2,NULL,XRES);
-		else
-			do_plot6(xx,y3,y1,NULL,y0,y2,NULL,XRES);
+		else {
+			if (outs == icSigCmykData)
+				do_plot6(xx,y3,y1,NULL,y0,y2,NULL,XRES);
+			else 	/* Assume RGB */
+				do_plot6(xx,NULL,y0,y1,y2,NULL,NULL,XRES);
+		}
 	}
 
 	/* Done with lookup object */
diff --git a/xicc/xcal.c b/xicc/xcal.c
index 06d343c..4400745 100644
--- a/xicc/xcal.c
+++ b/xicc/xcal.c
@@ -301,6 +301,9 @@ static int xcal_read(xcal *p, char *filename) {
 		return p->errc;
 	}
 
+	if (tcg->ntables < 1)
+		return 1;
+
 	rv = xcal_read_cgats(p, tcg, table, filename);
 
 	tcg->del(tcg);
diff --git a/xicc/xcam.c b/xicc/xcam.c
index 8e1a0cd..6117cd1 100644
--- a/xicc/xcam.c
+++ b/xicc/xcam.c
@@ -126,7 +126,7 @@ double Yb,		/* Relative Luminance of Background to reference white */
 double Lv,		/* Luminance of white in the Viewing/Scene/Image field (cd/m^2) */
 				/* Ignored if Ev is set to other than vc_none */
 double Yf,		/* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
-double Yg,		/* Glare as a fraction of the ambient (Y range 0.0 .. 1.0) */
+double Yg,		/* Glare as a fraction of the adapting/surround (Y range 0.0 .. 1.0) */
 double Gxyz[3],	/* The Glare white coordinates (typically the Ambient color) */
 int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
 ) {
diff --git a/xicc/xcam.h b/xicc/xcam.h
index 7ec6949..021c621 100644
--- a/xicc/xcam.h
+++ b/xicc/xcam.h
@@ -49,7 +49,7 @@ struct _icxcam {
 		double Lv,		/* Luminance of white in the Viewing/Scene/Image field (cd/m^2) */
 						/* Ignored if Ev is set */
 		double Yf,		/* Flare as a fraction of the reference white (range 0.0 .. 1.0) */
-		double Yg,		/* Glare as a fraction of the ambient (range 0.0 .. 1.0) */
+		double Yg,		/* Glare as a fraction of the adapting/surround (range 0.0 .. 1.0) */
 		double Gxyz[3],	/* The Glare white coordinates (typically the Ambient color) */
 		int hk			/* Flag, NZ to use Helmholtz-Kohlraush effect */
 	);
diff --git a/xicc/xfbview.c b/xicc/xfbview.c
index 6abf176..affdb53 100644
--- a/xicc/xfbview.c
+++ b/xicc/xfbview.c
@@ -25,6 +25,8 @@
 #include "numlib.h"
 #include "icc.h"
 #include "xicc.h"
+#include "vrml.h"
+#include "ui.h"
 
 #define RW 0.5		/* Device Delta */ 
 
@@ -65,7 +67,7 @@ void usage(void) {
 	fprintf(stderr," -d          Show PCS target -> average of device ref clippped PCS\n");
 	fprintf(stderr," -b          Show PCS target -> B2A lookup clipped PCS\n");
 	fprintf(stderr," -e          Show reference cliped PCS -> B2A lookup clipped PCS\n");
-	fprintf(stderr," -r res      Resolution of test grid\n");
+	fprintf(stderr," -r res      Resolution of test grid [Def 33]\n");
 	fprintf(stderr," -g          Do full grid, not just L = 0\n");
 	fprintf(stderr," -c          Do all values, not just clipped ones\n");
 	fprintf(stderr," -l tlimit   set total ink limit, 0 - 400%% (estimate by default)\n");
@@ -181,7 +183,7 @@ main(
 	strcpy(out_name, in_name);
 	if ((xl = strrchr(out_name, '.')) == NULL)	/* Figure where extention is */
 		xl = out_name + strlen(out_name);
-	strcpy(xl,".wrl");
+	xl[0] = '\000';				/* Remove extension */
 
 	/* Open up the file for reading */
 	if ((rd_fp = new_icmFileStd_name(in_name,"r")) == NULL)
@@ -201,7 +203,7 @@ main(
 		xicc *xicco;
 		icxLuBase *luo;
 		icxInk ink;					/* Ink parameters */
-		FILE *wrl;
+		vrml *wrl;
 		struct {
 			double x, y, z;
 			double wx, wy, wz;
@@ -250,60 +252,21 @@ main(
 			error("Expecting CMYK device");
 		}
 		
-		if ((wrl = fopen(out_name,"w")) == NULL) {
-			fprintf(stderr,"Error opening output file '%s'\n",out_name);
+		if ((wrl = new_vrml(out_name, doaxes, vrml_lab)) == NULL) {
+			fprintf(stderr,"new_vrml failed for '%s%s'\n",out_name,vrml_ext());
 			return 2;
 		}
 
-		/* Spit out a VRML 2 Object surface of gamut */
-
-		fprintf(wrl,"#VRML V2.0 utf8\n");
-		fprintf(wrl,"\n");
-		fprintf(wrl,"# Created by the Argyll CMS\n");
-		fprintf(wrl,"Transform {\n");
-		fprintf(wrl,"children [\n");
-		fprintf(wrl,"	NavigationInfo {\n");
-		fprintf(wrl,"		type \"EXAMINE\"        # It's an object we examine\n");
-		fprintf(wrl,"	} # We'll add our own light\n");
-		fprintf(wrl,"\n");
-		fprintf(wrl,"    DirectionalLight {\n");
-		fprintf(wrl,"        direction 0 0 -1      # Light illuminating the scene\n");
-		fprintf(wrl,"        direction 0 -1 0      # Light illuminating the scene\n");
-		fprintf(wrl,"    }\n");
-		fprintf(wrl,"\n");
-		fprintf(wrl,"    Viewpoint {\n");
-		fprintf(wrl,"        position 0 0 340      # Position we view from\n");
-		fprintf(wrl,"    }\n");
-		fprintf(wrl,"\n");
-		if (doaxes != 0) {
-			fprintf(wrl,"# Lab axes as boxes:\n");
-			for (i = 0; i < 5; i++) {
-				fprintf(wrl,"Transform { translation %f %f %f\n", axes[i].x, axes[i].y, axes[i].z);
-				fprintf(wrl,"\tchildren [\n");
-				fprintf(wrl,"\t\tShape{\n");
-				fprintf(wrl,"\t\t\tgeometry Box { size %f %f %f }\n",
-				                  axes[i].wx, axes[i].wy, axes[i].wz);
-				fprintf(wrl,"\t\t\tappearance Appearance { material Material ");
-				fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[i].r, axes[i].g, axes[i].b);
-				fprintf(wrl,"\t\t}\n");
-				fprintf(wrl,"\t]\n");
-				fprintf(wrl,"}\n");
-			}
-			fprintf(wrl,"\n");
-		}
-
 		/* ---------------------------------------------- */
 		/* The PCS target -> Reference clipped vectors */
 
 		if (doref) {
+			double rgb[3];
+
 			if (verb)
 				printf("Doing PCS target to reference clipped PCS Vectors\n");
 
-			fprintf(wrl,"\n");
-			fprintf(wrl,"Shape {\n");
-			fprintf(wrl,"  geometry IndexedLineSet { \n");
-			fprintf(wrl,"    coord Coordinate { \n");
-			fprintf(wrl,"	   point [\n");
+			wrl->start_line_set(wrl, 0);
 
 			i = 0;
 			for (coa[0] = 0; coa[0] < tres; coa[0]++) {
@@ -339,8 +302,8 @@ main(
 							printf("."), fflush(stdout);
 
 						/* Input PCS to ideal (Inverse AtoB) clipped PCS values */
-						fprintf(wrl,"%f %f %f,\n",in[1], in[2], in[0]-GAMUT_LCENT);
-						fprintf(wrl,"%f %f %f,\n",out[1], out[2], out[0]-GAMUT_LCENT);
+						wrl->add_vertex(wrl, 0, in);
+						wrl->add_vertex(wrl, 0, out);
 						i++;
 					}
 				}
@@ -350,17 +313,18 @@ main(
 
 			if (verb)
 				printf("\n");
-			fprintf(wrl,"      ]\n");
-			fprintf(wrl,"    }\n");
-			fprintf(wrl,"  coordIndex [\n");
 
 			for (j = 0; j < i; j++) {
-				fprintf(wrl,"%d, %d, -1,\n", j * 2, j * 2 + 1);
+				int ix[2];
+				ix[0] = j * 2;
+				ix[1] = j * 2 +1;
+				wrl->add_line(wrl, 0, ix);
 			}
-			fprintf(wrl,"    ]\n");
-			fprintf(wrl,"  }\n");
-			fprintf(wrl,"appearance Appearance { material Material { emissiveColor 1.0 0.1 0.1} }\n");
-			fprintf(wrl,"} # end shape\n");
+
+			rgb[0] = 1.0;
+			rgb[1] = 0.1;
+			rgb[2] = 0.1;
+			wrl->make_lines_cc(wrl, 0, 0.0, rgb);
 		}
 
 		/* ---------------------------------------------- */
@@ -368,14 +332,12 @@ main(
 		/* The PCS target -> clipped from average of surrounding device values, vectors */
 
 		if (dodelta) {
+			double rgb[3];
+
 			if (verb)
 				printf("Doing target PCS to average of 4 surrounding device to PCS Vectors\n");
 
-			fprintf(wrl,"\n");
-			fprintf(wrl,"Shape {\n");
-			fprintf(wrl,"  geometry IndexedLineSet { \n");
-			fprintf(wrl,"    coord Coordinate { \n");
-			fprintf(wrl,"	   point [\n");
+			wrl->start_line_set(wrl, 0);
 
 			i = 0;
 			for (coa[0] = 0; coa[0] < tres; coa[0]++) {
@@ -485,8 +447,8 @@ main(
 						if (verb)
 							printf("."), fflush(stdout);
 
-						fprintf(wrl,"%f %f %f,\n",in[1], in[2], in[0]-GAMUT_LCENT);
-						fprintf(wrl,"%f %f %f,\n",out[1], out[2], out[0]-GAMUT_LCENT);
+						wrl->add_vertex(wrl, 0, in);
+						wrl->add_vertex(wrl, 0, out);
 						i++;
 					}
 				}
@@ -496,26 +458,30 @@ main(
 
 			if (verb)
 				printf("\n");
-			fprintf(wrl,"      ]\n");
-			fprintf(wrl,"    }\n");
-			fprintf(wrl,"  coordIndex [\n");
 
 			for (j = 0; j < i; j++) {
-				fprintf(wrl,"%d, %d, -1,\n", j * 2, j * 2 + 1);
+				int ix[2];
+				ix[0] = j * 2;
+				ix[1] = j * 2 +1;
+				wrl->add_line(wrl, 0, ix);
 			}
-			fprintf(wrl,"    ]\n");
-			fprintf(wrl,"  }\n");
-			fprintf(wrl,"appearance Appearance { material Material { emissiveColor 0.9 0.9 0.9} }\n");
-			fprintf(wrl,"} # end shape\n");
 
+			rgb[0] = 0.9;
+			rgb[1] = 0.9;
+			rgb[2] = 0.9;
+			wrl->make_lines_cc(wrl, 0, 0.0, rgb);
 		}
 
 		/* ---------------------------------------------- */
 		/* The target PCS -> clipped PCS using B2A table vectore */
 
 		if (dob2a) {
+			double rgb[3];
+
 			icxLuBase *luoB;
 
+			wrl->start_line_set(wrl, 0);
+
 			/* Get a PCS to Device conversion object */
 			if ((luoB = xicco->get_luobj(xicco, ICX_CLIP_NEAREST, icmBwd, icAbsoluteColorimetric,
 			                             icSigLabData, icmLuOrdNorm, NULL, &ink)) == NULL) {
@@ -527,12 +493,6 @@ main(
 			if (verb)
 				printf("Doing target PCS to B2A clipped PCS Vectors\n");
 
-			fprintf(wrl,"\n");
-			fprintf(wrl,"Shape {\n");
-			fprintf(wrl,"  geometry IndexedLineSet { \n");
-			fprintf(wrl,"    coord Coordinate { \n");
-			fprintf(wrl,"	   point [\n");
-
 			i = 0;
 			for (coa[0] = 0; coa[0] < tres; coa[0]++) {
 				for (coa[1] = 0; coa[1] < tres; coa[1]++) {
@@ -563,8 +523,8 @@ main(
 						if (verb)
 							printf("."), fflush(stdout);
 
-						fprintf(wrl,"%f %f %f,\n",in[1], in[2], in[0]-GAMUT_LCENT);
-						fprintf(wrl,"%f %f %f,\n",out[1], out[2], out[0]-GAMUT_LCENT);
+						wrl->add_vertex(wrl, 0, in);
+						wrl->add_vertex(wrl, 0, out);
 						i++;
 					}
 				}
@@ -574,17 +534,18 @@ main(
 
 			if (verb)
 				printf("\n");
-			fprintf(wrl,"      ]\n");
-			fprintf(wrl,"    }\n");
-			fprintf(wrl,"  coordIndex [\n");
 
 			for (j = 0; j < i; j++) {
-				fprintf(wrl,"%d, %d, -1,\n", j * 2, j * 2 + 1);
+				int ix[2];
+				ix[0] = j * 2;
+				ix[1] = j * 2 +1;
+				wrl->add_line(wrl, 0, ix);
 			}
-			fprintf(wrl,"    ]\n");
-			fprintf(wrl,"  }\n");
-			fprintf(wrl,"appearance Appearance { material Material { emissiveColor 0.9 0.9 0.9} }\n");
-			fprintf(wrl,"} # end shape\n");
+
+			rgb[0] = 0.9;
+			rgb[1] = 0.9;
+			rgb[2] = 0.9;
+			wrl->make_lines_cc(wrl, 0, 0.0, rgb);
 
 			luoB->del(luoB);
 		}
@@ -593,8 +554,12 @@ main(
 		/* The reference clipped PCS -> B2A clipped PCS vectore */
 
 		if (doeee) {
+			double rgb[3];
+
 			icxLuBase *luoB;
 
+			wrl->start_line_set(wrl, 0);
+
 			/* Get a PCS to Device conversion object */
 			if ((luoB = xicco->get_luobj(xicco, ICX_CLIP_NEAREST, icmBwd, icAbsoluteColorimetric,
 			                             icSigLabData, icmLuOrdNorm, NULL, &ink)) == NULL) {
@@ -606,12 +571,6 @@ main(
 			if (verb)
 				printf("Doing reference clipped PCS to B2A table clipped PCS Vectors\n");
 
-			fprintf(wrl,"\n");
-			fprintf(wrl,"Shape {\n");
-			fprintf(wrl,"  geometry IndexedLineSet { \n");
-			fprintf(wrl,"    coord Coordinate { \n");
-			fprintf(wrl,"	   point [\n");
-
 			i = 0;
 			for (coa[0] = 0; coa[0] < tres; coa[0]++) {
 				for (coa[1] = 0; coa[1] < tres; coa[1]++) {
@@ -653,8 +612,8 @@ main(
 						if (verb)
 							printf("."), fflush(stdout);
 
-						fprintf(wrl,"%f %f %f,\n",check[1], check[2], check[0]-GAMUT_LCENT);
-						fprintf(wrl,"%f %f %f,\n",out[1], out[2], out[0]-GAMUT_LCENT);
+						wrl->add_vertex(wrl, 0, check);
+						wrl->add_vertex(wrl, 0, out);
 						i++;
 					}
 				}
@@ -664,31 +623,29 @@ main(
 
 			if (verb)
 				printf("\n");
-			fprintf(wrl,"      ]\n");
-			fprintf(wrl,"    }\n");
-			fprintf(wrl,"  coordIndex [\n");
 
 			for (j = 0; j < i; j++) {
-				fprintf(wrl,"%d, %d, -1,\n", j * 2, j * 2 + 1);
+				int ix[2];
+				ix[0] = j * 2;
+				ix[1] = j * 2 +1;
+				wrl->add_line(wrl, 0, ix);
 			}
-			fprintf(wrl,"    ]\n");
-			fprintf(wrl,"  }\n");
-			fprintf(wrl,"appearance Appearance { material Material { emissiveColor 0.9 0.9 0.9} }\n");
-			fprintf(wrl,"} # end shape\n");
+
+			rgb[0] = 0.9;
+			rgb[1] = 0.9;
+			rgb[2] = 0.9;
+			wrl->make_lines_cc(wrl, 0, 0.0, rgb);
 
 			luoB->del(luoB);
 		}
 
 		/* ---------------------------------------------- */
 
-		fprintf(wrl,"\n");
-		fprintf(wrl,"  ] # end of children for world\n");
-		fprintf(wrl,"}\n");
-	
-		if (fclose(wrl) != 0) {
-			fprintf(stderr,"Error closing output file '%s'\n",out_name);
+		if (wrl->flush(wrl) != 0) {
+			fprintf(stderr,"Error closing output file '%s%s'\n",out_name,vrml_ext());
 			return 2;
 		}
+		wrl->del(wrl);
 
 		/* Done with lookup object */
 		luo->del(luo);
diff --git a/xicc/xfit.c b/xicc/xfit.c
index d0912f9..d12919f 100644
--- a/xicc/xfit.c
+++ b/xicc/xfit.c
@@ -83,7 +83,8 @@
 							/* This seems to work badly, even with high smoothness. Why ? */
 							/* It does speed up 1D lut creation though. */
 
-#undef DEBUG 			/* Verbose debug information */
+#undef DEBUG 			/* Debug information */
+#undef DEBUG_PROGRESS 	/* Show powell progress */
 #undef DEBUG_PLOT 		/* Plot in & out curves */
 #undef SPECIAL_FORCE	/* Check rspl nodes against linear XYZ model */
 #undef SPECIAL_FORCE_GAMMA		/* Force correct gamma shaper curves */
@@ -99,8 +100,10 @@
 #define CURVEPOW 1.0    	/* Power to raise deltaE squared to in setting in/out curves */
 							/* This provides a means of punishing high maximum errors. */
 
-#define POWTOL 1e-4			/* Shaper Powell optimiser tollerance in delta E squared ^ CURVEPOW */
-#define MAXITS 2000			/* Shaper number of itterations before giving up */
+#define POWTOL1 1e-3		/* Shaper Powell optimiser tollerance for first passes */
+#define MAXITS1 1000		/* Shaper number of itterations for first passes */
+#define POWTOL 1e-5			/* Shaper Powell optimiser tollerance in delta E squared ^ CURVEPOW */
+#define MAXITS 4000			/* Shaper number of itterations before giving up */
 #define PDDEL  1e-6			/* Fake partial derivative del */
 
 /* Weights for shaper in/out curve parameters, to minimise unconstrained "wiggles" */
@@ -654,6 +657,27 @@ static void xfit_abs_to_rel(xfit *p, double *out, double *in) {
 	}
 }
 
+/* Convert an XYZ output value from absolute */
+/* to cLut relative using the current white point. */
+static void xfit_XYZ_abs_to_rel(xfit *p, double *out, double *in) {
+	if (p->flags & XFIT_OUT_WP_REL) {
+		if (p->flags & XFIT_OUT_LAB) {
+			icmMulBy3x3(out, p->fromAbs, in);
+			icmXYZ2Lab(&icmD50, out, out);
+		} else {
+			icmMulBy3x3(out, p->fromAbs, in);
+		}
+	} else {
+		if (p->flags & XFIT_OUT_LAB) {
+			icmXYZ2Lab(&icmD50, out, in);
+		} else {
+			out[0] = in[0];
+			out[1] = in[1];
+			out[2] = in[2];
+		}
+	}
+}
+
 /* - - - - - - - - - */
 
 /* return a weighting for the magnitude of the in and out */
@@ -921,7 +945,7 @@ static double xfitfunc(void *edata, double *v) {
 		smv = pow(smv, CURVEPOW);
 	rv = ev + smv;
 
-#ifdef DEBUG
+#ifdef DEBUG_PROGRESS
 if (xfitfunc_trace)
 fprintf(stdout,"~1(sm %f, ev %f)xfitfunc returning %f\n",smv,ev,rv);
 #endif
@@ -968,7 +992,6 @@ static double dxfitfunc(void *edata, double *dv, double *v) {
 
 	} else {
 		for (i = 0; i < p->opt_cnt; i++) { 
-//printf("~1 param %d = %f\n",i,v[i]);
 			p->v[p->opt_off + i] = v[i];
 		}
 	}
@@ -1112,13 +1135,12 @@ static double dxfitfunc(void *edata, double *dv, double *v) {
 	rv = ev + smv;
 
 	/* Sum the del for parameters being optimised and copy to return array */
-
 	if (p->opt_ssch) {
 		for (i = 0; i < p->sm_iluord; i++)
 			dv[i] = 0.0;
 		for (e = 0; e < di; e++) {	/* Combine per channel curve de's */
 			for (i = 0; i < p->sm_iluord; i++)
-				dv[i] += dav[p->shp_offs[e] + i] = sdav[p->shp_offs[e] + i];
+				dv[i] += dav[p->shp_offs[e] + i] + sdav[p->shp_offs[e] + i];
 		}
 		for (i = p->sm_iluord; i < p->opt_cnt; i++)	/* matrix and rest de's */ 
 			dv[i] = dav[p->mat_off + i - p->sm_iluord] + sdav[p->mat_off + i - p->sm_iluord];
@@ -1128,7 +1150,7 @@ static double dxfitfunc(void *edata, double *dv, double *v) {
 			dv[i] = dav[p->opt_off + i] + sdav[p->opt_off + i];
 	}
 
-#ifdef DEBUG
+#ifdef DEBUG_PROGRESS
 fprintf(stdout,"~1(sm %f, ev %f)dxfitfunc returning %f\n",smv,ev,rv);
 #endif
 
@@ -1251,7 +1273,7 @@ static double symoptfunc(void *edata, double *v) {
 
 	rv = out[0] * out[0];
 
-#ifdef DEBUG
+#ifdef DEBUG_PROGRESS
 printf("~1symoptfunc returning %f\n",rv);
 #endif
 	return rv;
@@ -1639,7 +1661,7 @@ printf("~1 changing %f %f %f -> %f %f %f\n", out[0], out[1], out[2], tout[0], to
 /* Do the fitting. */
 /* return nz on error */
 /* 1 = malloc or other error */
-int xfit_fit(
+static int xfit_fit(
 	struct _xfit *p, 
 	int flags,				/* Flag values */
 	int di,					/* Input dimensions */
@@ -1659,6 +1681,7 @@ int xfit_fit(
 	int gres[MXDI],			/* clut resolutions being optimised for/returned */
 	double out_min[MXDO],	/* Output value scaling/range minimum */
 	double out_max[MXDO],	/* Output value scaling/range maximum */
+//	co *bpo,				/* If != NULL, black point override in same spaces as ipoints */
 	double smooth,			/* clut rspl smoothing factor */
 	double oavgdev[MXDO],	/* Average output value deviation */
 	double demph,			/* dark emphasis factor for cLUT grid res. */
@@ -1678,6 +1701,9 @@ int xfit_fit(
 	double *b;				/* Base of parameters for this section */
 	int poff;
 
+	double powtol = POWTOL1;		/* powell/conjgrad initial tollerance */
+	int maxits = MAXITS1;			/* powell/conjgrad initial maximum itterations */
+
 	if (tcomb & oc_io)		/* If we're doing anything, we need the matrix */
 		tcomb |= oc_m;
 
@@ -1791,10 +1817,16 @@ int xfit_fit(
 		icmAry2XYZ(_wp, p->wp);
 
 		/* Absolute->Aprox. Relative Adaptation matrix */
-		icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
+		if (p->picc != NULL)
+			p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, icmD50, _wp, p->fromAbs);
+		else
+			icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
 	
 		/* Aproximate relative to absolute conversion matrix */
-		icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
+		if (p->picc != NULL)
+			p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, _wp, icmD50, p->toAbs);
+		else
+			icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
 
 		if (p->verb) {
 			double lab[3];
@@ -1921,11 +1953,11 @@ dump_xfit(p);
 		setup_xfit(p, p->wv, p->sa, 0.0, 0.5); 
 
 #ifdef NODDV
-		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 		                                 xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
 			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
 #else
-		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 		                       xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
 			warning("xfit_fit: Conjgrad failed to converge, residual error = %f", rerr);
 #endif
@@ -1935,14 +1967,58 @@ dump_xfit(p);
 #ifdef DEBUG
 printf("\nAfter matrix opt:\n");
 dump_xfit(p);
+
 #endif
 	}
 
 	/* Optimise input and matrix together */
 	if ((p->tcomb & oc_im) == oc_im) {
 		double rerr;
+		int sm_iluord = p->sm_iluord;
 
 		if (p->verb)
+			printf("About to optimise a common ord 0 input curve and matrix\n");
+
+		/* Setup pseudo-inverse if we need it */
+		if (p->flags & XFIT_FM_INPUT)
+			setup_piv(p);
+
+		p->opt_ssch = 1;
+		p->sm_iluord = 1;		/* Do a single order for first up */
+		p->opt_ch = -1;
+		p->opt_msk = oc_im;
+		setup_xfit(p, p->wv, p->sa, 0.5, 0.3); 
+
+#ifdef NODDV
+		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
+		                                xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
+#ifdef DEBUG
+			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
+#endif
+		}
+#else
+		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
+		                     xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
+#ifdef DEBUG
+			warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
+#endif
+		}
+#endif	/* !NODDV */
+		for (e = 0; e < di; e++) {				/* Copy optimised values back */
+			for (i = 0; i < p->sm_iluord; i++)
+				p->v[p->shp_offs[e] + i] = p->wv[i];
+			for (; i < p->iluord[e]; i++)
+				p->v[p->shp_offs[e] + i] = 0.0;
+		}
+		for (i = p->sm_iluord; i < p->opt_cnt; i++) 
+			p->v[p->mat_off + i - p->sm_iluord] = p->wv[i];
+#ifdef DEBUG
+printf("\nAfter single input and matrix opt:\n");
+dump_xfit(p);
+#endif
+
+		/* - - - - - - - - - - - */
+		if (p->verb)
 			printf("About to optimise a common input curve and matrix\n");
 
 		/* Setup pseudo-inverse if we need it */
@@ -1950,19 +2026,20 @@ dump_xfit(p);
 			setup_piv(p);
 
 		p->opt_ssch = 1;
+		p->sm_iluord = sm_iluord;		/* restore this */
 		p->opt_ch = -1;
 		p->opt_msk = oc_im;
 		setup_xfit(p, p->wv, p->sa, 0.5, 0.3); 
 
 #ifdef NODDV
-		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 		                                xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
 #ifdef DEBUG
 			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
 #endif
 		}
 #else
-		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 		                     xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
 #ifdef DEBUG
 			warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
@@ -1978,7 +2055,7 @@ dump_xfit(p);
 		for (i = p->sm_iluord; i < p->opt_cnt; i++) 
 			p->v[p->mat_off + i - p->sm_iluord] = p->wv[i];
 #ifdef DEBUG
-printf("\nAfter input and matrix opt:\n");
+printf("\nAfter single input and matrix opt:\n");
 dump_xfit(p);
 #endif
 
@@ -1986,6 +2063,11 @@ dump_xfit(p);
 		if (p->verb)
 			printf("About to optimise input curves and matrix\n");
 
+		if ((p->tcomb & oc_mo) != oc_mo) {	/* If this will be last fit */
+			powtol = POWTOL;
+			maxits = MAXITS;
+		}
+
 		/* Setup pseudo-inverse if we need it */
 		if (p->flags & XFIT_FM_INPUT)
 			setup_piv(p);
@@ -1998,14 +2080,14 @@ dump_xfit(p);
 		/* itterations and move on to the output curve, and worry about it not */
 		/* converging the second time through. */
 #ifdef NODDV
-		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 		                                xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
 #ifdef DEBUG
 			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
 #endif
 		}
 #else
-		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 		                     xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
 #ifdef DEBUG
 			warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
@@ -2027,6 +2109,11 @@ dump_xfit(p);
 		if (p->verb)
 			printf("About to optimise output curves and matrix\n");
 
+		if ((p->tcomb & oc_im) != oc_im) {	/* If this will be last fit */
+			powtol = POWTOL;
+			maxits = MAXITS;
+		}
+
 		/* Setup pseudo-inverse if we need it */
 		if (p->flags & XFIT_FM_INPUT)
 			setup_piv(p);
@@ -2036,11 +2123,11 @@ dump_xfit(p);
 		p->opt_msk = oc_mo;
 		setup_xfit(p, p->wv, p->sa, 0.3, 0.3); 
 #ifdef NODDV
-		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+		if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 		                                    xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
 			warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
 #else
-		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS, xfitfunc,
+		if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits, xfitfunc,
 		                                   dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
 			warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
 #endif
@@ -2057,6 +2144,15 @@ dump_xfit(p);
 			if (p->verb)
 				printf("About to optimise input curves and matrix again\n");
 
+
+#ifndef NODDV
+			if ((p->tcomb & oc_imo) != oc_imo)	/* If this will be last fit */
+#endif
+			{
+				powtol = POWTOL;
+				maxits = MAXITS;
+			}
+
 			/* Setup pseudo-inverse if we need it */
 			if (p->flags & XFIT_FM_INPUT)
 				setup_piv(p);
@@ -2066,11 +2162,11 @@ dump_xfit(p);
 			p->opt_msk = oc_im;
 			setup_xfit(p, p->wv, p->sa, 0.2, 0.2); 
 #ifdef NODDV
-			if (powell(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+			if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 			                               xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
 				warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
 #else
-			if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
+			if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
 			                           xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
 				warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
 #endif
@@ -2084,6 +2180,7 @@ dump_xfit(p);
 
 #ifndef NODDV
 		/* Optimise all together */
+		/* (This is very slow using powell) */
 		if ((p->tcomb & oc_imo) == oc_imo) {
 
 			if (p->verb)
@@ -2566,7 +2663,10 @@ printf("~1 ipos[%d][%d] = %f\n",e,i,cv);
 
 			/* Matrix needed to correct approx rel wp to target D50 */
 			icmAry2XYZ(_wp, wcc.v);		/* Aprox relative target white point */
-			icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->cmat);	/* Correction */
+			if (p->picc != NULL)	/* Correction */
+				p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, icmD50, _wp, p->cmat);
+			else
+				icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->cmat);
 
 			/* Compute the actual white point, and return it to caller */
 			icmMulBy3x3(wp, p->toAbs, wcc.v);
@@ -2586,8 +2686,13 @@ printf("~1 ipos[%d][%d] = %f\n",e,i,cv);
 
 			/* Fix absolute conversions to leave absolute response unchanged. */
 			icmAry2XYZ(_wp, wp);		/* Actual white point */
-			icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
-			icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
+			if (p->picc != NULL) {
+				p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, icmD50, _wp, p->fromAbs);
+				p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, _wp, icmD50, p->toAbs);
+			} else {
+				icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
+				icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
+			}
 
 			if (p->verb) {
 				double labwp[3];
@@ -2732,6 +2837,22 @@ printf("~1 ipos[%d][%d] = %f\n",e,i,cv);
 			);
 		}
 
+		/* Force black point to given value */
+//		if (bpo != NULL) {
+//			co tv;
+//			int rv;
+//
+//			xfit_inpscurves(p, tv.p, bpo->p);
+//
+//			xfit_XYZ_abs_to_rel(p, tv.v, bpo->v);
+//			xfit_invoutcurves(p, tv.v, tv.v);
+//printf("~1 xfit: fine after curves  black at %f %f %f to %f %f %f\n",
+//tv.p[0], tv.p[1], tv.p[2], tv.v[0], tv.v[1], tv.v[2]);
+//			rv = p->clut->tune_value(p->clut, &tv);
+//			if (rv != 0)
+//				warning("Black Point Override failed - clipping");
+//		}
+
 #ifdef SPECIAL_FORCE
 		/* Replace the rspl nodes with ones directly computed */
 		/* from the synthetic linear RGB->XYZ model */
@@ -2966,6 +3087,7 @@ static void xfit_del(xfit *p) {
 /* Create a transform fitting object */
 /* return NULL on error */
 xfit *new_xfit(
+icc *picc			/* ICC profile used to set cone space matrix, NULL for Bradford. */
 ) {
 	xfit *p;
 
@@ -2973,6 +3095,8 @@ xfit *new_xfit(
 		return NULL;
 	}
 
+	p->picc = picc;
+
 	/* Set method pointers */
 	p->fit         = xfit_fit;
 	p->incurve     = xfit_inpscurve;
diff --git a/xicc/xfit.h b/xicc/xfit.h
index 519f071..14b4ecc 100644
--- a/xicc/xfit.h
+++ b/xicc/xfit.h
@@ -70,6 +70,7 @@ typedef struct {
 
 /* Context for optimising input and output luts */
 struct _xfit {
+	icc *picc;				/* ICC profile used to set cone space matrix, NULL for Bradford. */
 	int verb;				/* Verbose */
 	int flags;				/* Behaviour flags */
 	int di, fdi;			/* Dimensionaluty of input and output */
@@ -170,6 +171,7 @@ struct _xfit {
 		int gres[MXDI],			/* clut resolutions being optimised for/returned */
 		double out_min[MXDO],	/* Output value scaling/range minimum */
 		double out_max[MXDO],	/* Output value scaling/range maximum */
+//		co *bpo,				/* If != NULL, black point override in same spaces as ipoints */
 		double smooth,			/* clut rspl smoothing factor */
 		double oavgdev[MXDO],	/* Average output value deviation */
 		double demph,			/* dark emphasis factor for cLUT grid res. */
@@ -200,7 +202,8 @@ struct _xfit {
 
 }; typedef struct _xfit xfit;
 
-xfit *new_xfit();
+/* The icc is to provide the cone space matrix. If NULL, Bradford will be used. */
+xfit *new_xfit(icc *picc);
 
 #endif /* XFIT_H */
 
diff --git a/xicc/xicc.c b/xicc/xicc.c
index a1c4531..a7556d5 100644
--- a/xicc/xicc.c
+++ b/xicc/xicc.c
@@ -6,7 +6,7 @@
  * Date:    2/7/00
  * Version: 1.00
  *
- * Copyright 2000, 2001 Graeme W. Gill
+ * Copyright 2000, 2001, 2014 Graeme W. Gill
  * All rights reserved.
  * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
  * see the License.txt file for licencing details.
@@ -60,7 +60,9 @@ icxLuBase * xicc_get_luobj(xicc *p, int flags, icmLookupFunc func, icRenderingIn
 static icxLuBase *xicc_set_luobj(xicc *p, icmLookupFunc func, icRenderingIntent intent,
                             icmLookupOrder order, int flags, int no, int nobw, cow *points,
 							icxMatrixModel *skm,
-                            double dispLuminance, double wpscale, double smooth, double avgdev,
+                            double dispLuminance, double wpscale,
+//							double *bpo,
+                            double smooth, double avgdev,
 							double demph, icxViewCond *vc, icxInk *ink, xcal *cal, int quality);
 static void icxLutSpaces(icxLuBase *p, icColorSpaceSignature *ins, int *inn,
                          icColorSpaceSignature *outs, int *outn,
@@ -680,12 +682,6 @@ double *kblack			/* XYZ Output. Looked up if possible or set to black[] otherwis
 #ifdef DEBUG
 		printf("~1 Lab pivot %f %f %f, Lab K direction %f %f %f\n",bfs.p1[0],bfs.p1[1],bfs.p1[2],bfs.p2[0],bfs.p2[1],bfs.p2[2]);
 #endif
-		/* Start with the K only as the current best value */
-		brv = bpfindfunc((void *)&bfs, dblack);
-#ifdef DEBUG
-		printf("~1 initial brv for K only = %f\n",brv);
-#endif
-
 		/* Set the random start 0 location as 000K */
 		/* and the random start 1 location as CMY0 */
 		{
@@ -707,6 +703,12 @@ double *kblack			/* XYZ Output. Looked up if possible or set to black[] otherwis
 			rs1[kch] = 0.0;		/* K value */
 		}
 
+		/* Start with the K only as the current best value */
+		brv = bpfindfunc((void *)&bfs, dblack);
+#ifdef DEBUG
+		printf("~1 initial brv for K only = %f\n",brv);
+#endif
+
 		/* Find the device black point using optimization */
 		/* Do several trials to avoid local minima. */
 		rand32(0x12345678);	/* Make trial values deterministic */
@@ -1056,6 +1058,7 @@ cow *points,				/* Array of input points in target PCS space */
 icxMatrixModel *skm,   		/* Optional skeleton model (used for input profiles) */
 double dispLuminance,		/* > 0.0 if display luminance value and is known */
 double wpscale,				/* > 0.0 if input white point is to be scaled */
+//double *bpo,				/* != NULL for black point override XYZ */
 double smooth,				/* RSPL smoothing factor, -ve if raw */
 double avgdev,				/* reading Average Deviation as a proportion of the input range */
 double demph,				/* dark emphasis factor for cLUT grid res. */
@@ -1104,12 +1107,17 @@ int quality					/* Quality metric, 0..3 */
     	case icmMatrixFwdType:
 			if (smooth < 0.0)
 				smooth = -smooth;
-			xplu = set_icxLuMatrix(p, plu, flags, no, nobw, points, skm, dispLuminance, wpscale, quality, smooth);
+			xplu = set_icxLuMatrix(p, plu, flags, no, nobw, points, skm, dispLuminance, wpscale,
+//			                       bpo,
+			                       quality, smooth);
 			break;
 
     	case icmLutType:
 			/* ~~~ Should add check that it is a fwd profile ~~~ */
-			xplu = set_icxLuLut(p, plu, func, intent, flags, no, nobw, points, skm, dispLuminance, wpscale, smooth, avgdev, demph, vc, ink, quality);
+			xplu = set_icxLuLut(p, plu, func, intent, flags, no, nobw, points, skm, dispLuminance,
+			                    wpscale,
+//			                    bpo,
+			                    smooth, avgdev, demph, vc, ink, quality);
 			break;
 
 		default:
@@ -1213,14 +1221,14 @@ icxViewCond *vc		/* Viewing parameters to return */
 	/* Numbers we're trying to find */
 	ViewingCondition Ev = vc_none;
 	double Wxyz[3] = {-1.0, -1.0, -1.0};	/* Adapting white color */
-	double La = -1.0;		/* Adapting luminance */
+	double La = -1.0;						/* Adapting/Surround luminance */
 	double Ixyz[3] = {-1.0, -1.0, -1.0};	/* Illuminant color */
 	double Li = -1.0;						/* Illuminant luminance */
 	double Lb = -1.0;		/* Backgrount luminance */
 	double Yb = -1.0;		/* Background relative luminance to Lv */
 	double Lve = -1.0;		/* Emissive device image luminance */
 	double Lvr = -1.0;		/* Reflective device image luminance */
-	double Lv = -1.0;		/* device image luminance */
+	double Lv = -1.0;		/* Device image luminance */
 	double Yf = -1.0;		/* Flare relative luminance to Lv */
 	double Yg = -1.0;		/* Glare relative luminance to La */
 	double Gxyz[3] = {-1.0, -1.0, -1.0};	/* Glare color */
@@ -1909,7 +1917,7 @@ icxViewCond *vc
 	if (vc->Ev == vc_none)
 		printf("  Image luminance = %f cd/m^2\n",vc->Lv);
 	printf("  Flare to image ratio = %f\n",vc->Yf);
-	printf("  Glare to ambient ratio = %f\n",vc->Yg);
+	printf("  Glare to adapting/surround ratio = %f\n",vc->Yg);
 	printf("  Flare color = %f %f %f\n",vc->Gxyz[0], vc->Gxyz[1], vc->Gxyz[2]);
 }
 
@@ -2150,7 +2158,7 @@ char *as				/* Alias string selector, NULL for none */
 		gmi->bph = gmm_bendBP;		/* extent and bend */
 		gmi->gamcpf  = 1.0;			/* Full gamut compression */
 		gmi->gamexf  = 0.0;			/* No gamut expansion */
-		gmi->gamcknf  = 0.8;		/* High Sigma knee in gamut compress */
+		gmi->gamcknf  = 0.9;		/* 0.9 High Sigma knee in gamut compress */
 		gmi->gamxknf  = 0.0;		/* No knee in gamut expand */
 		gmi->gampwf  = 1.0;			/* Full Perceptual surface weighting factor */
 		gmi->gamswf  = 0.0;			/* No Saturation surface weighting factor */
@@ -2176,7 +2184,7 @@ char *as				/* Alias string selector, NULL for none */
 		gmi->bph = gmm_bendBP;		/* extent and bend */
 		gmi->gamcpf  = 1.0;			/* Full gamut compression */
 		gmi->gamexf  = 0.0;			/* No gamut expansion */
-		gmi->gamcknf  = 0.8;		/* High Sigma knee in gamut compress */
+		gmi->gamcknf  = 0.9;		/* 0.9 High Sigma knee in gamut compress */
 		gmi->gamxknf  = 0.0;		/* No knee in gamut expand */
 		gmi->gampwf  = 1.0;			/* Full Perceptual surface weighting factor */
 		gmi->gamswf  = 0.0;			/* No Saturation surface weighting factor */
@@ -3616,8 +3624,7 @@ double *in			/* Input di values */
 /* Including partial derivative for input and parameters. */
 
 
-/* 3x3 matrix multiplication, with the matrix in a 1D array */
-/* with respect to the input and parameters. */
+/* 3x3 matrix in 1D array multiplication */
 void icxMulBy3x3Parm(
 	double out[3],			/* Return input multiplied by matrix */
 	double mat[9],			/* Matrix organised in [slow][fast] order */
@@ -3639,7 +3646,39 @@ void icxMulBy3x3Parm(
 }
 
 
-/* 3x3 matrix multiplication, with partial derivatives */
+/* 3x3 matrix in 1D array multiplication, with partial derivatives */
+/* with respect to just the input. */
+void icxdpdiiMulBy3x3Parm(
+	double out[3],			/* Return input multiplied by matrix */
+	double din[3][3],		/* Return deriv for each [output] with respect to [input] */
+	double mat[9],			/* Matrix organised in [slow][fast] order */
+	double in[3]			/* Input values */
+) {
+	double *v, ov[3];
+	int e, f;
+
+	/* Compute the output values */
+	v = mat;
+	for (f = 0; f < 3; f++) {
+		ov[f] = 0.0;						/* For each output value */
+		for (e = 0; e < 3; e++) {
+			ov[f] += *v++ * in[e];
+		}
+	}
+
+	/* Compute deriv. with respect to the input values */
+	/* This is pretty simple for a matrix ... */
+	v = mat;
+	for (f = 0; f < 3; f++)
+		for (e = 0; e < 3; e++)
+			din[f][e] = *v++;
+
+	out[0] = ov[0];
+	out[1] = ov[1];
+	out[2] = ov[2];
+}
+
+/* 3x3 matrix in 1D array multiplication, with partial derivatives */
 /* with respect to the input and parameters. */
 void icxdpdiMulBy3x3Parm(
 	double out[3],			/* Return input multiplied by matrix */
@@ -3683,163 +3722,6 @@ void icxdpdiMulBy3x3Parm(
 	out[2] = ov[2];
 }
 
-/* ------------------------------------------- */
-/* BT.1886 support */
-
-/* Compute technical gamma from effective gamma in BT.1886 style */
-
-/* Info for optimization */
-typedef struct {
-	double thyr;		/* 50% input target */
-	double roo;			/* 0% input target */
-} gam_fits;
-
-/* gamma + input offset function handed to powell() */
-static double gam_fit(void *dd, double *v) {
-	gam_fits *gf = (gam_fits *)dd;
-	double gamma = v[0];
-	double a, b;
-	double rv = 0.0;
-	double tt;
-
-	if (gamma < 0.0) {
-		rv += 100.0 * -gamma;
-		gamma = 1e-4;
-	}
-
-	tt = pow(gf->roo, 1.0/gamma);
-	b = tt/(1.0 - tt);						/* Offset */
-	a = pow(1.0 - tt, gamma);				/* Gain */
-
-	tt = a * pow((0.5 + b), gamma);
-	tt = tt - gf->thyr;
-	rv += tt * tt;
-	
-	return rv;
-}
-
-/* Given the effective gamma and the output offset Y, */
-/* return the technical gamma needed for the correct 50% response. */
-double xicc_tech_gamma(
-	double egamma,			/* effective gamma needed */
-	double off				/* Output offset required */
-) {
-	gam_fits gf;
-	double op[1], sa[1], rv;
-
-	if (off <= 0.0) {
-		return egamma;
-	}
-
-	gf.thyr = pow(0.5, egamma);					/* Advetised 50% target */
-	gf.roo = off;
-
-	op[0] = egamma;
-	sa[0] = 0.1;
-
-	if (powell(&rv, 1, op, sa, 1e-6, 500, gam_fit, (void *)&gf, NULL, NULL) != 0)
-		warning("Computing effective gamma and input offset is inaccurate");
-
-	return op[0];
-}
-
-
-/* Set the bt1886_info to a default do nothing state */
-void bt1886_setnop(bt1886_info *p) {
-	p->ingo = 0.0;
-	p->outsc = 1.0;
-	p->outL = 0.0;
-	p->tab[0] = 0.0;
-	p->tab[1] = 0.0;
-}
-
-/* Setup the bt1886_info for the given target */
-void bt1886_setup(bt1886_info *p, double *XYZbp, double gamma) {
-	double Lab[3], bkipow;
-	p->gamma = gamma;
-
-	icmXYZ2Lab(&icmD50, Lab, XYZbp);
-
-	p->outL = Lab[0];		/* For bp blend */
-	p->tab[0] = Lab[1];		/* a* b* correction needed */
-	p->tab[1] = Lab[2];
-
-	bkipow = pow(XYZbp[1], 1.0/p->gamma);
-	p->ingo = bkipow/(1.0 - bkipow);		/* non-linear Y that makes out black point */
-	p->outsc = pow(1.0 - bkipow, p->gamma);	/* Scale to restore 1 -> 1 */
-}
-
-/* Apply BT.1886 black offset and gamma curve to the XYZ out of the input profile. */
-/* Do this in the colorspace defined by the input profile matrix lookup, */
-/* so it will be relative XYZ. We assume that BT.1886 does a Rec709 to gamma */
-/* viewing adjustment, irrespective of the source profile transfer curve. */
-void bt1886_apply(bt1886_info *p, icmLuMatrix *lu, double *out, double *in) {
-	int j;
-	double vv;
-
-#ifdef DEBUG
-	printf("bt1886 XYZ in %f %f %f\n", in[0],in[1],in[2]);
-#endif
-
-	lu->bwd_matrix(lu, out, in);
-
-#ifdef DEBUG
-	printf("bt1886 RGB in %f %f %f\n", out[0],out[1],out[2]);
-#endif
-
-	for (j = 0; j < 3; j++) {
-		vv = out[j];
-	
-		/* Convert linear light to Rec709 transfer curve */
-		if (vv < 0.018)
-			vv = 4.5 * vv;
-		else
-			vv = 1.099 * pow(vv, 0.45) - 0.099;
-		
-		/* Apply input offset & re-scale, and then gamma of 2.4/custom gamma */
-		vv = vv + p->ingo;
-		
-		if (vv > 0.0)
-			vv = p->outsc * pow(vv, p->gamma);
-
-		out[j] = vv;
-	}
-
-	lu->fwd_matrix(lu, out, out);
-
-#ifdef DEBUG
-	printf("bt1886 RGB bt.1886 %f %f %f\n", out[0],out[1],out[2]);
-#endif
-
-	icmXYZ2Lab(&icmD50, out, out);
-
-#ifdef DEBUG
-	printf("bt1886 Lab after Y adj. %f %f %f\n", out[0],out[1],out[2]);
-#endif
-
-	/* Blend ab to required black point offset p->tab[] as L approaches black. */
-	vv = (out[0] - p->outL)/(100.0 - p->outL);	/* 0 at bp, 1 at wp */
-	vv = 1.0 - vv;
-
-	if (vv < 0.0)
-		vv = 0.0;
-	else if (vv > 1.0)
-		vv = 1.0;
-	vv = pow(vv, 40.0);
-	out[1] += vv * p->tab[0];
-	out[2] += vv * p->tab[1];
-
-#ifdef DEBUG
-	printf("bt1886 Lab after wp adj. %f %f %f\n", out[0],out[1],out[2]);
-#endif
-
-	icmLab2XYZ(&icmD50, out, out);
-
-#ifdef DEBUG
-	printf("bt1886 XYZ out %f %f %f\n", out[0],out[1],out[2]);
-#endif
-}
-
 /* - - - - - - - - - - */
 #undef stricmp
 
@@ -3863,5 +3745,3 @@ void bt1886_apply(bt1886_info *p, icmLuMatrix *lu, double *out, double *in) {
 
 
 
-
-
diff --git a/xicc/xicc.h b/xicc/xicc.h
index 2e69ef1..37fff5d 100644
--- a/xicc/xicc.h
+++ b/xicc/xicc.h
@@ -101,6 +101,7 @@ const char *icx2str(icmEnumType etype, int enumval);
 struct _icxMatrixModel {
 	void *imp;		/* Opaque implementation */
 
+	icc *picc;		/* ICC profile used to set cone space matrix, NULL for Bradford. */
 	int isLab;		/* Convert lookup to Lab */
 
 	void (*force) (struct _icxMatrixModel *p, double *targ, double *in);
@@ -112,6 +113,7 @@ struct _icxMatrixModel {
 /* Create a matrix model of a set of points, and return an object to lookup */
 /* points from the model. Return NULL on error. */
 icxMatrixModel *new_MatrixModel(
+icc *picc,			/* ICC profile used to set cone space matrix, NULL for Bradford. */
 int verb,			/* NZ if verbose */
 int nodp,			/* Number of points */
 cow *ipoints,		/* Array of input points in XYZ space */
@@ -186,7 +188,7 @@ typedef struct {
 	double Lv;			/* Luminance of white in the Image/Scene/Viewing field (cd/m^2) */
 						/* Ignored if Ev is set to other than vc_none */
 	double Yf;			/* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
-	double Yg;			/* Glare as a fraction of the ambient (Y range 0.0 .. 1.0) */
+	double Yg;			/* Glare as a fraction of the adapting/surround (Y range 0.0 .. 1.0) */
 	double Gxyz[3];		/* The Glare white coordinates (ie the Ambient color) */
 						/* will be taken from Wxyz if Gxyz <= 0.0 */
 	char *desc;			/* Possible description of this VC */
@@ -302,6 +304,7 @@ struct _xicc {
 	                                                                    /* value and is known */
 	                                  double wpscale,			/* > 0.0 if input white pt is */
 	                                                                    /* is to be scaled */
+//									  double *bpo, 				/* != NULL black point override */
 	                                  double smooth,			/* RSPL smoothing factor, */
 	                                                                        /* -ve if raw */
 	                                  double avgdev,			/* Avge Dev. of points */
@@ -897,15 +900,23 @@ double *in			/* Input di values */
 
 /* - - - - - - - - - - */
 
-/* 3x3 matrix multiplication, with the matrix in a 1D array */
-/* with respect to the input and parameters. */
+/* 3x3 matrix in 1D array multiplication */
 void icxMulBy3x3Parm(
 	double out[3],			/* Return input multiplied by matrix */
 	double mat[9],			/* Matrix organised in [slow][fast] order */
 	double in[3]			/* Input values */
 );
 
-/* 3x3 matrix multiplication, with partial derivatives */
+/* 3x3 matrix in 1D array multiplication, with partial derivatives */
+/* with respect to just the input. */
+void icxdpdiiMulBy3x3Parm(
+	double out[3],			/* Return input multiplied by matrix */
+	double din[3][3],		/* Return deriv for each [output] with respect to [input] */
+	double mat[9],			/* Matrix organised in [slow][fast] order */
+	double in[3]			/* Input values */
+);
+
+/* 3x3 matrix in 1D array multiplication, with partial derivatives */
 /* with respect to the input and parameters. */
 void icxdpdiMulBy3x3Parm(
 	double out[3],			/* Return input multiplied by matrix */
@@ -915,32 +926,6 @@ void icxdpdiMulBy3x3Parm(
 	double in[3]			/* Input values */
 );
 
-/* ------------------------------------------- */
-/* BT.1886 support */
-
-/* Convert an effective gamma given an offset into a technical gamma */
-double xicc_tech_gamma(double egamma, double off);
-
-typedef struct {
-	double ingo;				/* input Y gamma offset for bt1886 */
-	double outsc;				/* output Y scale for bt1886 */
-	double outL;				/* output black point L value */
-	double tab[2];				/* Target ab offset value at zero input for bt1886 */
-	double gamma;				/* bt.1886 technical gamma to apply */
-} bt1886_info;
-
-/* Set the bt1886_info to a default do nothing state */
-void bt1886_setnop(bt1886_info *p);
-
-/* Setup the bt1886_info for the given target */
-void bt1886_setup(bt1886_info *p, double *XYZbp, double gamma);
-
-/* Apply BT.1886 black offset and gamma curve to */
-/* the XYZ out of the input profile. */
-/* Do this in the colorspace defined by the input profile matrix lookup, */
-/* so it will be relative XYZ */
-void bt1886_apply(bt1886_info *p, icmLuMatrix *lu, double *out, double *in);
-
 /* - - - - - - - - - - */
 
 #include "xcal.h"
diff --git a/xicc/xicclu.c b/xicc/xicclu.c
index 4987e65..f3b1267 100644
--- a/xicc/xicclu.c
+++ b/xicc/xicclu.c
@@ -34,8 +34,9 @@
 #include "copyright.h"
 #include "aconfig.h"
 #include "numlib.h"
-#include "plot.h"
 #include "xicc.h"
+#include "plot.h"
+#include "ui.h"
 
 #undef SPTEST				/* Test rspl gamut surface code */
 
@@ -152,7 +153,7 @@ void spioutf(void *cbntx, double *out, double *in) {
 
 int
 main(int argc, char *argv[]) {
-	int fa,nfa;				/* argument we're looking at */
+	int fa, nfa, mfa;				/* argument we're looking at */
 	char prof_name[MAXNAMEL+1];
 	icmFile *fp = NULL;
 	icc *icco = NULL;
@@ -187,8 +188,8 @@ main(int argc, char *argv[]) {
 	int repLCh = 0;			/* Report LCh */
 	int repXYZ100 = 0;		/* Scale XYZ by 10 */
 	double scale = 0.0;		/* Device value scale factor */
-	int in_tvenc;			/* 1 to use RGB Video Level encoding, 2 = Rec601, 3 = Rec709 YCbCr */
-	int out_tvenc;			/* 1 to use RGB Video Level encoding, 2 = Rec601, 3 = Rec709 YCbCr */
+	int in_tvenc = 0;		/* 1 to use RGB Video Level encoding, 2 = Rec601, 3 = Rec709 YCbCr */
+	int out_tvenc = 0;		/* 1 to use RGB Video Level encoding, 2 = Rec601, 3 = Rec709 YCbCr */
 	int rv = 0;
 	char buf[200];
 	double uin[MAX_CHAN], in[MAX_CHAN], out[MAX_CHAN], uout[MAX_CHAN];
@@ -224,7 +225,8 @@ main(int argc, char *argv[]) {
 		usage("Too few arguments");
 
 	/* Process the arguments */
-	for(fa = 1;fa < argc;fa++) {
+	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 */
 			char *na = NULL;		/* next argument after flag, null if none */
@@ -232,7 +234,7 @@ main(int argc, char *argv[]) {
 			if (argv[fa][2] != '\000')
 				na = &argv[fa][2];		/* next is directly after flag */
 			else {
-				if ((fa+1) < argc) {
+				if ((fa+1+mfa) < argc) {
 					if (argv[fa+1][0] != '-') {
 						nfa = fa + 1;
 						na = argv[nfa];		/* next is seperate non-flag argument */
@@ -244,11 +246,10 @@ main(int argc, char *argv[]) {
 				usage("Requested usage");
 
 			/* Verbosity */
-			else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') {
-				fa = nfa;
-				if (na == NULL)
+			else if (argv[fa][1] == 'v') {
+				if (na == NULL) {
 					verb = 2;
-				else {
+				} else {
 					if (na[0] == '0')
 						verb = 0;
 					else if (na[0] == '1')
@@ -257,6 +258,7 @@ main(int argc, char *argv[]) {
 						verb = 2;
 					else
 						usage("Illegal verbosity level");
+					fa = nfa;
 				}
 			}
 
@@ -266,7 +268,6 @@ main(int argc, char *argv[]) {
 			}
 			/* Plot start or end override */
 			else if (argv[fa][1] == 'G') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -G");
 				if (na[0] == 's' || na[0] == 'S') {
 					if (sscanf(na+1,":%lf:%lf:%lf",&pstart[0],&pstart[1],&pstart[2]) != 3)
@@ -276,6 +277,7 @@ main(int argc, char *argv[]) {
 						usage("Unrecognised parameters after -Ge");
 				} else
 					usage("Unrecognised parameters after -G");
+				fa = nfa;
 			}
 			/* Actual target values */
 			else if (argv[fa][1] == 'a') {
@@ -299,10 +301,10 @@ main(int argc, char *argv[]) {
 			}
 			/* Device scale */
 			else if (argv[fa][1] == 's') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -s");
 				scale = atof(na);
 				if (scale <= 0.0) usage("Illegal scale value");
+				fa = nfa;
 			}
 			/* Video RGB encoding */
 			else if (argv[fa][1] == 'e'
@@ -343,7 +345,6 @@ main(int argc, char *argv[]) {
 
 			/* function */
 			else if (argv[fa][1] == 'f') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -f");
     			switch (na[0]) {
 					case 'f':
@@ -375,11 +376,11 @@ main(int argc, char *argv[]) {
 					default:
 						usage("Unknown parameter after flag -f");
 				}
+				fa = nfa;
 			}
 
 			/* Intent */
 			else if (argv[fa][1] == 'i') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -i");
     			switch (na[0]) {
 					case 'p':
@@ -405,11 +406,11 @@ main(int argc, char *argv[]) {
 					default:
 						usage("Unknown parameter after flag -i");
 				}
+				fa = nfa;
 			}
 
 			/* PCS override */
 			else if (argv[fa][1] == 'p') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -i");
     			switch (na[0]) {
 					case 'x':
@@ -465,11 +466,11 @@ main(int argc, char *argv[]) {
 					default:
 						usage("Unknown parameter after flag -i");
 				}
+				fa = nfa;
 			}
 
 			/* Search order */
 			else if (argv[fa][1] == 'o') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -o");
     			switch (na[0]) {
 					case 'n':
@@ -483,12 +484,12 @@ main(int argc, char *argv[]) {
 					default:
 						usage("Unknown parameter after flag -o");
 				}
+				fa = nfa;
 			}
 
 			/* Inking rule */
 			else if (argv[fa][1] == 'k'
 			      || argv[fa][1] == 'K') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -k");
 				if (argv[fa][1] == 'k')
 					locus = 0;			/* K value target */
@@ -573,18 +574,19 @@ main(int argc, char *argv[]) {
 					default:
 						usage("Unknown parameter after flag -k");
 				}
+				fa = nfa;
 			}
 
 			else if (argv[fa][1] == 'l') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -l");
 				tlimit = atoi(na)/100.0;
+				fa = nfa;
 			}
 
 			else if (argv[fa][1] == 'L') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -L");
 				klimit = atoi(na)/100.0;
+				fa = nfa;
 			}
 
 #ifdef SPTEST
@@ -597,7 +599,6 @@ main(int argc, char *argv[]) {
 #endif
 			/* Viewing conditions */
 			else if (argv[fa][1] == 'c') {
-				fa = nfa;
 				if (na == NULL) usage("No parameter after flag -c");
 #ifdef NEVER
 				if (na[0] >= '0' && na[0] <= '9') {
@@ -658,6 +659,7 @@ main(int argc, char *argv[]) {
 						usage("Unrecognised parameters after -cg");
 				} else
 					usage("Unrecognised parameters after -c");
+				fa = nfa;
 			}
 
 			else 
@@ -674,14 +676,6 @@ main(int argc, char *argv[]) {
 			error("chrom_locus_poligon failed");
 	}
 
-	if (verb > 1) {
-		icmFile *op;
-		if ((op = new_icmFileStd_fp(stdout)) == NULL)
-			error ("Can't open stdout");
-		icco->header->dump(icco->header, op, 1);
-		op->del(op);
-	}
-
 	/* Open up the profile for reading */
 	if ((fp = new_icmFileStd_name(prof_name,"r")) == NULL)
 		error ("Can't open file '%s'",prof_name);
@@ -959,6 +953,15 @@ main(int argc, char *argv[]) {
 		outn = inn = cal->devchan;
 	}
 
+	if (verb > 1 && icco != NULL) {
+		icmFile *op;
+		if ((op = new_icmFileStd_fp(stdout)) == NULL)
+			error ("Can't open stdout");
+		icco->header->dump(icco->header, op, 1);
+		op->del(op);
+	}
+
+
 	if (doplot) {
 		int i, j;
 		double xx[XRES];
diff --git a/xicc/xlut.c b/xicc/xlut.c
index bbad934..5b07ca5 100644
--- a/xicc/xlut.c
+++ b/xicc/xlut.c
@@ -6,7 +6,7 @@
  * Date:    2/7/00
  * Version: 1.00
  *
- * Copyright 2000, 2001 Graeme W. Gill
+ * Copyright 2000, 2001, 2014 Graeme W. Gill
  * All rights reserved.
  * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
  * see the License.txt file for licencing details.
@@ -155,6 +155,8 @@
 #undef REPORT_LOCUS_SEGMENTS    /* [Undef[ Examine how many segments there are in aux inversion */
 
 #define XYZ_EXTRA_SMOOTH 20.0		/* Extra smoothing factor for XYZ profiles */
+									/* !!! Note this is mainly due to smoothing being */
+									/* scaled by data range in rspl code !!! */
 #define SHP_SMOOTH 1.0	/* Input shaper curve smoothing */
 #define OUT_SMOOTH1 1.0	/* Output shaper curve smoothing for L*, X,Y,Z */
 #define OUT_SMOOTH2 1.0	/* Output shaper curve smoothing for a*, b* */
@@ -177,13 +179,6 @@
  *       error(), should return status.
  */
 
-#ifndef _CAT2
-#define _CAT2(n1,n2)  n1 ## n2
-#define CAT2(n1,n2) _CAT2(n1,n2)
-#endif
-
-
-
 static double icxLimitD(icxLuLut *p, double *in);		/* For input' */
 #define icxLimitD_void ((double (*)(void *, double *))icxLimitD)	/* Cast with void 1st arg */
 static double icxLimit(icxLuLut *p, double *in);		/* For input */
@@ -2855,6 +2850,7 @@ cow                *ipoints,		/* Array of input points (Lab or XYZ normalized to
 icxMatrixModel     *skm,    		/* Optional skeleton model (used for input profiles) */
 double             dispLuminance,	/* > 0.0 if display luminance value and is known */
 double             wpscale,			/* > 0.0 if white point is to be scaled */
+//double            *bpo,				/* != NULL for XYZ black point override dev & XYZ */
 double             smooth,			/* RSPL smoothing factor, -ve if raw */
 double             avgdev,			/* reading Average Deviation as a prop. of the input range */
 double             demph,			/* dark emphasis factor for cLUT grid res. */
@@ -2878,6 +2874,7 @@ int                quality			/* Quality metric, 0..3 */
 	double oavgdev[MXDO];	/* Average output value deviation */
 	int gres[MXDI];			/* RSPL/CLUT resolution */
 	xfit *xf = NULL;		/* Curve fitting class instance */
+//	co bpop;				/* bpo dev + XYZ value */
 
 	if (flags & ICX_VERBOSE)
 		rsplflags |= RSPL_VERBOSE;
@@ -2933,6 +2930,9 @@ int                quality			/* Quality metric, 0..3 */
 	/* very "straight", and the lack of tension reduces any noise reduction effect. */
 	/* !!! This probably means that we should switch to 3rd order smoothness criteria !! */
 	/* We apply an arbitrary correction here */
+	/* !!!! There is also a bug in the rspl code, where smoothness is */
+	/* scaled by data range. This is making Lab smoothing ~100 times */
+	/* more than XYZ smoothing. Fix this with SMOOTH2 changes ?? */
 	if (p->pcs == icSigXYZData) {
 		oavgdev[0] = XYZ_EXTRA_SMOOTH * 0.70 * avgdev;
 		oavgdev[1] = XYZ_EXTRA_SMOOTH * 1.00 * avgdev;
@@ -3176,6 +3176,8 @@ int                quality			/* Quality metric, 0..3 */
 							nw++;
 						}
 					}
+					/* Setup bpo device value in case we need it */
+//					bpop.p[0] = bpop.p[1] = bpop.p[2] = 0.0;
 					break;
 	
 				case icSigGrayData: {	/* Could be additive or subtractive */
@@ -3207,6 +3209,7 @@ int                quality			/* Quality metric, 0..3 */
 						nw = nminwp;
 						if (minwp[pcsy]/nminwp < (0.5 * pcsymax))
 							nw = 0;					/* Looks like a mistake */
+//						bpop.p[0] = 1.0;
 					}
 					if (nmaxwp > 0				/* Additive */
 					 && (nminwp == 0 || maxwp[pcsy]/nmaxwp > minwp[pcsy]/nminwp)) {
@@ -3216,6 +3219,7 @@ int                quality			/* Quality metric, 0..3 */
 						nw = nmaxwp;
 						if (maxwp[pcsy]/nmaxwp < (0.5 * pcsymax))
 							nw = 0;					/* Looks like a mistake */
+//						bpop.p[0] = 0.0;
 					}
 					break;
 				}
@@ -3240,6 +3244,12 @@ int                quality			/* Quality metric, 0..3 */
 			wp[2] /= (double)nw;
 			if (p->pcs != icSigXYZData) 	/* Convert white point to XYZ */
 				icmLab2XYZ(&icmD50, wp, wp);
+
+//			if (bpo != NULL) {			/* Copy black override XYZ value */
+//				bpop.v[0] = bpo[0];
+//				bpop.v[1] = bpo[1];
+//				bpop.v[2] = bpo[2];
+//			}
 		}
 
 		if (flags & ICX_VERBOSE) {
@@ -3274,7 +3284,7 @@ int                quality			/* Quality metric, 0..3 */
 
 		optcomb tcomb = oc_ipo;	/* Create all by default */
 
-		if ((xf = CAT2(new_, xfit)()) == NULL) {
+		if ((xf = new_xfit(icco)) == NULL) {
 			p->pp->errc = 2;
 			sprintf(p->pp->err,"Creation of xfit object failed");
 			p->del((icxLuBase *)p);
@@ -3374,6 +3384,7 @@ int                quality			/* Quality metric, 0..3 */
 		if (xf->fit(xf, xfflags, p->inputChan, p->outputChan,
 			rsplflags, wp, dwhite, wpscale, dgwhite, 
 		    ipoints, nodp, skm, in_min, in_max, gres, out_min, out_max,
+//			bpo != NULL ? &bpop : NULL,
 		    smooth, oavgdev, demph, iord, sord, oord, shp_smooth, out_smooth, tcomb,
 		   (void *)p, xfit_to_de2, xfit_to_dde2) != 0) {
 			p->pp->errc = 2;
@@ -3533,7 +3544,7 @@ int                quality			/* Quality metric, 0..3 */
 				/* to use for the rich black. */
 				for (e = 0; e < p->inputChan; e++)
 					bcc.p[e] = 0.0;
-				if (p->ink.klimit < 0.0)
+				if (p->ink.klimit <= 0.0)
 					bcc.p[kch] = 1.0;
 				else
 					bcc.p[kch] = p->ink.klimit;		/* K value */
@@ -3552,10 +3563,6 @@ int                quality			/* Quality metric, 0..3 */
 				if (flags & ICX_VERBOSE)
 					printf("K only black direction (Lab) =                      %f %f %f\n",bfs.p2[0], bfs.p2[1], bfs.p2[2]);
 #endif
-				/* Start with the K only as the current best value */
-				brv = bfindfunc((void *)&bfs, bcc.p);
-//printf("~1 initial brv for K only = %f\n",brv);
-
 				/* Set the random start 0 location as 000K */
 				/* and the random start 1 location as CMY0 */
 				{
@@ -3563,12 +3570,12 @@ int                quality			/* Quality metric, 0..3 */
 
 					for (e = 0; e < p->inputChan; e++)
 						rs0[e] = 0.0;
-					if (p->ink.klimit < 0.0)
+					if (p->ink.klimit <= 0.0)
 						rs0[kch] = 1.0;
 					else
 						rs0[kch] = p->ink.klimit;		/* K value */
 
-					if (p->ink.tlimit < 0.0)
+					if (p->ink.tlimit <= 0.0)
 						tv = 1.0;
 					else
 						tv = p->ink.tlimit/(p->inputChan - 1.0);
@@ -3577,6 +3584,12 @@ int                quality			/* Quality metric, 0..3 */
 					rs1[kch] = 0.0;		/* K value */
 				}
 
+				/* Start with the K only as the current best value */
+				for (e = 0; e < p->inputChan; e++)
+					bcc.p[e] = rs0[e];
+				brv = bfindfunc((void *)&bfs, bcc.p);
+//printf("~1 initial brv for K only = %f\n",brv);
+
 				/* Find the device black point using optimization */
 				/* Do several trials to avoid local minima. */
 				rand32(0x12345678);	/* Make trial values deterministic */
@@ -3697,6 +3710,13 @@ int                quality			/* Quality metric, 0..3 */
 				wp[i] *= scale;
 				bp[i] *= scale;
 			}
+
+//			if (bpo != NULL) {
+//				bp[0] = bpo[0];
+//				bp[1] = bpo[1];
+//				bp[2] = bpo[2];
+//				printf("Overide Black point XYZ = %s, Lab = %s\n", icmPdv(3,bp),icmPLab(bp));
+//			}
 		}
 
 		if (h->deviceClass == icSigDisplayClass
diff --git a/xicc/xmatrix.c b/xicc/xmatrix.c
index 53db237..034a396 100644
--- a/xicc/xmatrix.c
+++ b/xicc/xmatrix.c
@@ -34,6 +34,8 @@
  *
  */
 
+
+
 #define USE_CIE94_DE	/* Use CIE94 delta E measure when creating fit */
 
 /* Weights in shaper parameters, to minimise unconstrained "wiggles" */
@@ -51,6 +53,7 @@
 
 #undef DEBUG			/* Extra printfs */
 #undef DEBUG_PLOT		/* Plot curves */
+#define G_DEB 0			/* g_deb default value */ 
 
 /* ========================================================= */
 /* Forward and Backward Matrix type conversion */
@@ -84,9 +87,8 @@ double *in			/* Vector of input values */
 	int rv = 0;
 	rv |= ((icmLuMatrix *)p->plu)->fwd_abs((icmLuMatrix *)p->plu, out, in);
 
-	if (p->pcs == icxSigJabData) {
+	if (p->pcs == icxSigJabData)
 		p->cam->XYZ_to_cam(p->cam, out, out);
-	}
 	return rv;
 }
 
@@ -614,6 +616,8 @@ double *v			/* Pointer to parameters */
 	return XSHAPE_MAG * tparam/3.0;
 }
 
+int g_deb = G_DEB;
+
 /* Matrix optimisation function handed to powell() */
 static double mxoptfunc(void *edata, double *v) {
 	mxopt *p = (mxopt *)edata;
@@ -621,6 +625,8 @@ static double mxoptfunc(void *edata, double *v) {
 	double xyz[3], lab[3];
 	int i;
 
+	if (g_deb) printf("\n");
+
 	for (i = 0; i < p->nodp; i++) {
 
 		/* Apply our function */
@@ -629,7 +635,7 @@ static double mxoptfunc(void *edata, double *v) {
 	
 		/* Convert to Lab */
 		icmXYZ2Lab(&p->wp, lab, xyz);
-//printf("%f %f %f -> %f %f %f, target %f %f %f\n", p->points[i].p[0], p->points[i].p[1], p->points[i].p[2], lab[0], lab[1], lab[2], p->points[i].v[0], p->points[i].v[1], p->points[i].v[2]);
+if (g_deb) printf("%d: %f %f %f -> %f %f %f, target %f %f %f, w %f\n", i, p->points[i].p[0], p->points[i].p[1], p->points[i].p[2], lab[0], lab[1], lab[2], p->points[i].v[0], p->points[i].v[1], p->points[i].v[2],p->points[i].w);
 	
 		/* Accumulate total delta E squared */
 #ifdef USE_CIE94_DE
@@ -674,6 +680,7 @@ static double mxoptfunc(void *edata, double *v) {
 	rv += err * 1000.0;
 
 #ifdef DEBUG
+if (g_deb)
 printf("~9(%f)mxoptfunc returning %f\n",smv,rv);
 #endif
 
@@ -774,8 +781,8 @@ double scale		/* Scale device values */
 	/* Set quality/effort  factors */
 	if (quality >= 3) {			/* Ultra high */
 		os->norders = 20;
-		maxits = 10000;
-		stopon = 5e-7;
+		maxits = 50000;
+		stopon = 1e-14;
 	} else if (quality == 2) {	/* High */
 		os->norders = 12;
 		maxits = 5000;
@@ -845,15 +852,22 @@ double scale		/* Scale device values */
 				icmLab2XYZ(&icmD50, points[i].v, points[i].v);
 			icmXYZ2Lab(&os->wp, points[i].v, points[i].v);
 			icmLab2LCh(lch, points[i].v);
+
 			/* Apply any neutral weighting */
 			if (lch[1] < 10.0) {
 				double w = nweight;
 				if (lch[1] > 5.0)
 					w = 1.0 + (nweight - 1.0) * (10.0 - lch[1])/(10.0 - 5.0);
-				points[i].w = w;
+				points[i].w *= w;
 			}
 //printf("~1 patch %d = Lab %f %f %f, C = %f w = %f\n",i,points[i].v[0], points[i].v[1], points[i].v[2], lch[1],points[i].w);
 		}
+
+#if !defined(NOT_PRIVATE) && defined(HACK)
+# pragma message("!!!!!!!!!!!!!!! xicc/xmatrix.c HACK code enabled !!!!!!!!!!!!!!!!!!")
+		printf("!!!! HACK: setting white point ixt %d weight to zero\n",wix);
+		points[wix].w = 0.0;
+#endif
 	}
 
 	/* Set initial matrix optimisation values */
@@ -1025,6 +1039,7 @@ double scale		/* Scale device values */
 				if (os->verb)
 					printf("Creating matrix and curves...\n"); 
 		
+//g_deb = 1;
 				if (powell(&rerr, os->optdim, os->v, os->sa, stopon, maxits,
 				           mxoptfunc, (void *)os, mxprogfunc, (void *)os) != 0)
 					warning("Powell failed to converge, residual error = %f",rerr);
@@ -1117,7 +1132,10 @@ static void icxMM_force_exact(icxMatrixModel *p, double *targ, double *rgb) {
 
 	icmAry2XYZ(_ap, axyz);
 	icmAry2XYZ(_tp, txyz);
-	icmChromAdaptMatrix(ICM_CAM_BRADFORD, _tp, _ap, cmat);
+	if (p->picc != NULL)
+		p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, _tp, _ap, cmat);
+	else
+		icmChromAdaptMatrix(ICM_CAM_BRADFORD, _tp, _ap, cmat);
 
 	/* Apply correction to fine tune matrix. */
 	mxtransform(os, cmat);
@@ -1140,6 +1158,7 @@ static void icxMM_del(icxMatrixModel *p) {
 /* Create a matrix model of a set of points, and return an object to lookup */
 /* points from the model. Return NULL on error. */
 icxMatrixModel *new_MatrixModel(
+icc *picc,			/* ICC profile used to set cone space matrix, NULL for Bradford. */
 int verb,			/* NZ if verbose */
 int nodp,			/* Number of points */
 cow *ipoints,		/* Array of input points in XYZ space */
@@ -1159,6 +1178,7 @@ double scale		/* Scale device values */
 	if ((p = (icxMatrixModel *) calloc(1,sizeof(icxMatrixModel))) == NULL)
 		return NULL;
 
+	p->picc = picc;
 	p->force = icxMM_force_exact;
 	p->lookup = icxMM_lookup;
 	p->del = icxMM_del;
@@ -1201,6 +1221,7 @@ cow                *ipoints,		/* Array of input points in XYZ space */
 icxMatrixModel     *skm,    		/* Optional skeleton model (not used here) */
 double             dispLuminance,	/* > 0.0 if display luminance value and is known */
 double             wpscale,			/* > 0.0 if input white point is to be scaled */
+//double             *bpo,			/* != NULL for XYZ black point override dev & XYZ */
 int                quality,			/* Quality metric, 0..3 */
 double             smooth			/* Curve smoothing, nominally 1.0 */
 ) {
@@ -1485,10 +1506,10 @@ double             smooth			/* Curve smoothing, nominally 1.0 */
 			icmAry2XYZ(_wp, wp);
 	
 			/* Absolute->Aprox. Relative Adaptation matrix */
-			icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, fromAbs);
+			icco->chromAdaptMatrix(icco, ICM_CAM_NONE, icmD50, _wp, fromAbs);
 		
 			/* Aproximate relative to absolute conversion matrix */
-			icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, toAbs);
+			icco->chromAdaptMatrix(icco, ICM_CAM_NONE, _wp, icmD50, toAbs);
 		}
 
 	} else {
@@ -1497,7 +1518,8 @@ double             smooth			/* Curve smoothing, nominally 1.0 */
 	}
 
 	/* Create copy of input points with output converted to white relative */
-	if ((rpoints = (cow *)malloc(nodp * sizeof(cow))) == NULL) {
+	/* Allow one extra point for possible bpo value */
+	if ((rpoints = (cow *)malloc((nodp+1) * sizeof(cow))) == NULL) {
 		xicp->errc = 1;
 		sprintf(xicp->err,"set_icxLuMatrix: malloc failed");
 		p->del((icxLuBase *)p);
@@ -1513,7 +1535,43 @@ double             smooth			/* Curve smoothing, nominally 1.0 */
 		/* abs out -> aprox. rel out */
 		icmMulBy3x3(rpoints[i].v, fromAbs, rpoints[i].v);
 	}
-  
+
+#ifdef NEVER
+	/* If black point override and shaper curves */
+	if (bpo != NULL && !isLinear && !isGamma) {
+		double tw = 0.0;		/* Total weight */
+
+printf("Got bpo\n");
+		/* Zero out any black data points, and sum up total weihting */
+		for (i = 0; i < nodp; i++) {
+			if (rpoints[i].p[0] < 0.001				/* We're assuming RGB */
+			 && rpoints[i].p[1] < 0.001
+			 && rpoints[i].p[2] < 0.001) {
+				rpoints[i].w = 0.0;
+printf("Zero'd point %d\n",i);
+			}
+			tw += rpoints[i].w;
+		}
+printf("Total weight = %f\n",tw);
+
+		/* Add our override black point */
+		/* and give it a dominant weighting */
+		for (e = 0; e < inputChan; e++)
+			rpoints[nodp].p[e] = 0.0;
+		for (f = 0; f < outputChan; f++)
+			rpoints[nodp].v[f] = bpo[f];
+printf(" set black to %f %f %f\n", bpo[0], bpo[1], bpo[2]);
+
+		/* abs out -> aprox. rel out */
+		icmMulBy3x3(rpoints[nodp].v, fromAbs, rpoints[nodp].v);
+
+		rpoints[nodp].w = 20.0 * tw;
+printf(" set black %d w = %f\n", nodp,rpoints[nodp].w);
+
+		nodp++;
+	}
+#endif // NEVER
+
 	/* ------------------------------- */
 
 	/* (Use a gamma curve as 0th order shape) */
@@ -1529,6 +1587,10 @@ double             smooth			/* Curve smoothing, nominally 1.0 */
 	}
 	free(rpoints); rpoints = NULL;
 
+#if !defined(NOT_PRIVATE) && defined(HACK)
+# pragma message("!!!!!!!!!!!!!!! xicc/xmatrix.c HACK code enabled !!!!!!!!!!!!!!!!!!")
+		printf("!!!! HACK: skipping white point fine tune\n");
+#else
 	/* The overall device to absolute conversion is now what we want */
 	/* (as dictated by the points, weighting and best fit), */
 	/* but we need to adjust the device to relative conversion */
@@ -1553,7 +1615,7 @@ double             smooth			/* Curve smoothing, nominally 1.0 */
 	
 		/* Matrix needed to correct aprox white to target D50 */
 		icmAry2XYZ(_wp, aw);		/* Aprox relative target white point */
-		icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, cmat);	/* Correction */
+		icco->chromAdaptMatrix(icco, ICM_CAM_NONE, icmD50, _wp, cmat);	/* Correction */
 	
 		/* Compute the current absolute white point */
 		icmMulBy3x3(wp, toAbs, aw);
@@ -1563,8 +1625,8 @@ double             smooth			/* Curve smoothing, nominally 1.0 */
 	
 		/* Fix relative conversions to leave absolute response unchanged. */
 		icmAry2XYZ(_wp, wp);		/* Actual white point */
-		icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, fromAbs);
-		icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, toAbs);
+		icco->chromAdaptMatrix(icco, ICM_CAM_NONE, icmD50, _wp, fromAbs);
+		icco->chromAdaptMatrix(icco, ICM_CAM_NONE, _wp, icmD50, toAbs);
 
 		if (flags & ICX_VERBOSE) {
 			double tw[3];
@@ -1573,6 +1635,7 @@ double             smooth			/* Curve smoothing, nominally 1.0 */
 			printf("                 abs WP = XYZ %s, Lab %s\n", icmPdv(3, wp), icmPLab(wp));
 		}
 	}
+#endif
 
 	/* Create default wpscale */
 	if (wpscale < 0.0) {
@@ -1657,8 +1720,8 @@ double             smooth			/* Curve smoothing, nominally 1.0 */
 
 		/* Fix absolute conversions to leave absolute response unchanged. */
 		icmAry2XYZ(_wp, wp);		/* Actual white point */
-		icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, fromAbs);
-		icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, toAbs);
+		icco->chromAdaptMatrix(icco, ICM_CAM_NONE, icmD50, _wp, fromAbs);
+		icco->chromAdaptMatrix(icco, ICM_CAM_NONE, _wp, icmD50, toAbs);
 	}
 
 	/* Look up the actual black point */
diff --git a/xicc/xspect.c b/xicc/xspect.c
index cc0ce85..477892b 100644
--- a/xicc/xspect.c
+++ b/xicc/xspect.c
@@ -99,6 +99,73 @@ static int gcc_bug_fix(int i) {
 /* Dummy "no illuminant" illuminant spectra used to signal an emmission */
 /* or equal energy 'E' illuminant */
 static xspect il_none = {
+	531, 300.0, 830.0,	/* 531 bands from 300 to 830 in 1nm steps */
+	1.0,				/* Scale factor */
+	{
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
+		1.0
+	}
+};
+ 
+#ifdef NEVER
+static xspect il_none = {
 	54, 300.0, 830.0,	/* 54 bands from 300 to 830 in 10nm steps */
 	1.0,				/* Scale factor */
 	{
@@ -110,6 +177,7 @@ static xspect il_none = {
 		1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
 	}
 };
+#endif /* NEVER */
  
 
 /* CIE 15.2-1986 Table 1.1 */
@@ -478,9 +546,13 @@ double temp					/* Optional temperature in degrees kelvin, for Dtemp and Ptemp *
 				uv_filter(&il_D50M2, &il_D50);
 			*sp = il_D50M2;
 			return 0;
+	    case icxIT_D55:
+			return daylight_il(sp, 5500.0);
 	    case icxIT_D65:
 			*sp = il_D65;
 			return 0;
+	    case icxIT_D75:
+			return daylight_il(sp, 7500.0);
 	    case icxIT_E:
 			*sp = il_none;
 			return 0;
@@ -3427,6 +3499,7 @@ void xspect_denorm(xspect *sp) {
 }
 
 #ifndef SALONEINSTLIB
+
 /* Convert from one xspect type to another (targ type) */
 /* Linear or polinomial interpolation will be used as appropriate */
 /* (converted to targ norm too) */
@@ -3458,6 +3531,54 @@ void xspect2xspect(xspect *dst, xspect *targ, xspect *src) {
 }
 
 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+/* Plot up to 3 spectra */
+void xspect_plot(xspect *sp1, xspect *sp2, xspect *sp3) {
+	double xx[XSPECT_MAX_BANDS];
+	double y1[XSPECT_MAX_BANDS];
+	double y2[XSPECT_MAX_BANDS];
+	double y3[XSPECT_MAX_BANDS];
+	int j;
+	double wl, wlshort, wllong;
+
+	if (sp1 == NULL)
+		return;
+
+	wlshort = sp1->spec_wl_short;
+	wllong = sp1->spec_wl_long;
+
+	if (sp2 != NULL) {
+		if (sp2->spec_wl_short < wlshort)
+			wlshort = sp2->spec_wl_short;
+		if (sp2->spec_wl_long > wllong)
+			wllong = sp2->spec_wl_long;
+	}
+	
+	if (sp3 != NULL) {
+		if (sp3->spec_wl_short < wlshort)
+			wlshort = sp3->spec_wl_short;
+		if (sp3->spec_wl_long > wllong)
+			wllong = sp3->spec_wl_long;
+	}
+	
+	wlshort = floor(wlshort + 0.5);
+	wllong = floor(wllong + 0.5);
+
+	/* Compute at 1nm intervals over the whole range covered */
+	for (j = 0, wl = wlshort; j < XSPECT_MAX_BANDS && wl < wllong; j++, wl += 1.0) {
+#if defined(__APPLE__) && defined(__POWERPC__)
+		gcc_bug_fix(j);
+#endif
+		xx[j] = wl;
+		y1[j] = value_xspect(sp1, wl);
+		if (sp2 != NULL)
+			y2[j] = value_xspect(sp2, wl);
+		if (sp3 != NULL)
+			y3[j] = value_xspect(sp3, wl);
+	}
+	do_plot(xx, y1, sp2 != NULL ? y2 : NULL, sp3 != NULL ? y3 : NULL, j);
+}
+
+/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
 /* Given an emission spectrum, set the UV output to the given level. */
 /* The shape of the UV is taken from FWA1_stim, and the level is */
@@ -3625,6 +3746,13 @@ static int xsp2cie_fwa_apply(xsp2cie *p, xspect *out, xspect *in);
 	FWA spectrum and the actual spectrum under the illuminant to
 	create the correction model. This could be fine tuned by doing
 	similar measurements of neutral patches.
+
+	Other possible limitations:
+
+	Instrument illuminant spectrum shape:
+	It is assumed it is stable and 'A' like. Aging of
+	the lamp may invalidate this assumption ?
+	
  */
 
 /*
@@ -3659,11 +3787,12 @@ static int xsp2cie_set_fwa_imp(xsp2cie *p) {
 
 	DBG("set_fwa started\n");
 
-	p->bw = 1.0;		/* Intergrate over 1nm bands */
+	p->bw = 1.0;				/* Intergrate over 1nm bands */
 	p->oillum = p->illuminant;	/* Take copy of observer illuminant */
 	xspect_denorm(&p->oillum);
-	if (p->tillum.spec_n == 0) {	/* If not set by set_fwa(), use observer illuminant */
+	if (p->tillum.spec_n == 0) {	/* If not set by set_fwa(), copy observer illuminant */
 		p->tillum = p->oillum;		/* as target/simulated instrument illuminant. */
+		DBG("using observer illum as FWA target\n");
 	}
 
 	/* Compute Y = 1 normalised instrument illuminant spectrum */
@@ -4026,7 +4155,7 @@ static int xsp2cie_set_fwa_imp(xsp2cie *p) {
 static int xsp2cie_set_fwa(xsp2cie *p,	/* this */
 xspect *iillum,		/* Spectrum of instrument illuminent */
 xspect *tillum,		/* Spectrum of target/simulated instrument illuminant */
-					/* NULL to use observer illuminant. */
+					/* NULL to use observer model illuminant. */
 xspect *media		/* Spectrum of plain media measured under that instrument */
 ) {
 	p->iillum = *iillum;			/* Take copy of instrument illuminant */
@@ -4035,26 +4164,32 @@ xspect *media		/* Spectrum of plain media measured under that instrument */
 		p->tillum = *tillum;		/* Take copy of target/simulated instrument illuminant */
 		xspect_denorm(&p->tillum);	/* Remove normalisation factor */
 	} else {
-		p->tillum.spec_n = 0;
+		p->tillum.spec_n = 0;		/* Use observer model illum. as FWA source */
 	}
 	p->imedia = *media;		/* Take copy of measured media */
 
 	return xsp2cie_set_fwa_imp(p);
 }
 
-/* Set FWA given updated conversion illuminant. */
+/* Set FWA given updated conversion illuminants. */
 /* We assume that xsp2cie_set_fwa has been called first. */
 static int xsp2cie_update_fwa_custillum(
 xsp2cie *p,	/* this */
 xspect *tillum,		/* Spectrum of target/simulated instrument illuminant, */
                     /* NULL to use previous set_fwa() value. */
-xspect *custIllum	/* Spectrum of observer illuminant */
+xspect *custIllum	/* Spectrum of observer model illuminant */
+					/* NULL to use previous new_xsp2cie() value. */
 ) {
 	if (tillum != NULL) {
 		p->tillum = *tillum;		/* Take copy of target/simulated instrument illuminant */
 		xspect_denorm(&p->tillum);	/* Remove normalisation factor */
+	} else {
+		p->tillum.spec_n = 0;		/* Use observer model illum. as FWA source */
+	}
+
+	if (custIllum != NULL) {
+		p->illuminant = *custIllum;		/* Updated observer model illuminant */
 	}
-	p->illuminant = *custIllum;
 
 	return xsp2cie_set_fwa_imp(p);
 }
@@ -4655,6 +4790,17 @@ void xsp2cie_convert(xsp2cie *p, double *out, xspect *in) {
 	xsp2cie_sconvert(p, NULL, out, in);
 }
 
+/* Return the illuminant XYZ being used in the CIE XYZ/Lab conversion. */ 
+/* Note that this will returne the 'E' illuminant XYZ for emissive. */
+void xsp2cie_get_cie_il(xsp2cie *p, double *xyz) {
+	xspect sp;
+
+	standardIlluminant(&sp, icxIT_E, 0.0); 
+	p->convert(p, xyz, &sp);
+	if (p->doLab)
+		icmLab2XYZ(&icmD50, xyz, xyz);
+}
+
 void xsp2cie_del(
 xsp2cie *p
 ) {
@@ -4684,7 +4830,7 @@ icxClamping clamp				/* NZ to clamp XYZ/Lab to be +ve */
 			p->isemis = 1;
 			break;
 	    case icxIT_custom:
-			p->illuminant = *custIllum;
+			p->illuminant = *custIllum;			/* Struct copy */
 			break;
 	    case icxIT_A:
 			p->illuminant = il_A;
@@ -4701,9 +4847,14 @@ icxClamping clamp				/* NZ to clamp XYZ/Lab to be +ve */
 				uv_filter(&il_D50M2, &il_D50);
 			p->illuminant = il_D50M2;
 			break;
+	    case icxIT_D55:
+			daylight_il(&p->illuminant, 5500.0);
+			break;
 	    case icxIT_D65:
 			p->illuminant = il_D65;
 			break;
+	    case icxIT_D75:
+			daylight_il(&p->illuminant, 7500.0);
 	    case icxIT_E:
 			p->illuminant = il_none;
 			break;
@@ -4789,6 +4940,7 @@ icxClamping clamp				/* NZ to clamp XYZ/Lab to be +ve */
 
 	p->convert      = xsp2cie_convert;
 	p->sconvert     = xsp2cie_sconvert;
+	p->get_cie_il   = xsp2cie_get_cie_il;
 #ifndef SALONEINSTLIB
 	p->set_mw       = xsp2cie_set_mw;		/* Default no media white */
 	p->set_fwa      = xsp2cie_set_fwa;		/* Default no FWA compensation */
@@ -4848,7 +5000,8 @@ int icx_spectrum_locus(double xyz[3], double wl, icxObserverType obType) {
 }
 
 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-/* Pre-calculated spectral locuses of Daylight and Plankian at 5 Mired intervals */
+/* Pre-calculated spectral locuses of Daylight and Plankian at 5 Mired intervals, */
+/* created using illlocus.c */
 /* These aren't actually spectrum, they are XYZ values */
 /* indexed by temperature in Mired */
 
@@ -5256,7 +5409,6 @@ typedef struct {
 	xspect *iloc;			/* Locus to match to */
 	double xyz[3];			/* Target XYZ */
 	icmXYZNumber XYZ;		/* Target as XYZ number for DE wp */
-	xsp2cie *conv;			/* Means of converting spectrum to XYZ */
 	int viscct;				/* nz to use visual best match color temperature */
 } cct2ctx;
 
@@ -5264,7 +5416,6 @@ static double cct2_func(void *fdata, double tp[]) {
 	cct2ctx *x = (cct2ctx *)fdata;
 	double xyz[3];		/* Current value */
 	double lab1[3], lab2[3];
-	xspect sp;
 	double rv = 0.0;
 	icmXYZNumber *wp = &x->XYZ;
 
@@ -5289,6 +5440,13 @@ static double cct2_func(void *fdata, double tp[]) {
 		rv = icmLabDEsq(lab1, lab2);
 	}
 	
+	/* Discourage going beyond ends of locus */
+	if (tp[0] < x->iloc->spec_wl_short ) {
+		rv += 5000.0 * (x->iloc->spec_wl_short - tp[0]);
+	} else if (tp[0] > x->iloc->spec_wl_long) {
+		rv += 5000.0 * (tp[0] - x->iloc->spec_wl_long);
+	} 
+
 //a1logd(g_log, 1, " cct2_func returning %f for temp = %f\n",rv,1e6/tp[0]);
 //DBGF((DBGA,"returning %f for temp = %f\n",rv,tp[0]));
 	return rv;
@@ -5362,7 +5520,6 @@ int viscct				/* nz to use visual CIEDE2000, 0 to use CCT CIE 1960 UCS. */
 
 	/* Locate the CCT in Mired */
 	if (powell(&rv, 1, cp, s, 0.01, 1000, cct2_func, (void *)&x, NULL, NULL) != 0) {
-		x.conv->del(x.conv);
 		return -1.0;
 	}
 
@@ -5379,6 +5536,60 @@ int viscct				/* nz to use visual CIEDE2000, 0 to use CCT CIE 1960 UCS. */
 	return 1e6/cp[0];
 }
 
+/* Given a choice of temperature dependent illuminant (icxIT_Dtemp or icxIT_Ptemp), */
+/* a color temperature and a Y value, return the corresponding XYZ */
+/* An observer type can be chosen for interpretting the spectrum of the input and */
+/* the illuminant. */
+/* Return xyz[0] = -1.0 on erorr */
+void icx_ill_ct2XYZ(
+double xyz[3],			/* Return the XYZ value */
+icxIllumeType ilType,	/* Type of illuminant, icxIT_Dtemp or icxIT_Ptemp */
+icxObserverType obType,	/* Observer, CIE_1931_2 or CIE_1964_10 */
+int viscct,				/* nz to use visual CIEDE2000, 0 to use CCT CIE 1960 UCS. */
+double tin,				/* Input temperature */
+double Yin				/* Input Y value */
+) {
+	xspect *iloc;			/* Locus to match to */
+
+	double cp[1], s[1];
+	
+	if (ilType != icxIT_Dtemp && ilType != icxIT_Ptemp) {
+		xyz[0] = -1.0;
+		return;
+	}
+	if (obType != icxOT_CIE_1931_2 && obType != icxOT_CIE_1964_10) {
+		xyz[0] = -1.0;
+		return;
+	}
+
+	/* Locus to use */
+	if (obType == icxOT_CIE_1931_2) {
+		if (ilType == icxIT_Dtemp) {
+			iloc = illoc_Daylight_CIE_1931_2;
+		} else {
+			iloc = illoc_Plankian_CIE_1931_2;
+		}
+	} else {
+		if (ilType == icxIT_Dtemp) {
+			iloc = illoc_Daylight_CIE_1964_10;
+		} else {
+			iloc = illoc_Plankian_CIE_1964_10;
+		}
+	}
+
+	/* Convert temperature to mired */
+	tin =  1e6/tin;
+
+	/* Get XYZ for given temp in Mired. */
+	/* Will clip to limits of locus */
+	getval_raw_xspec3_lin(iloc, xyz, tin);
+
+	/* Scale by Yin */
+	xyz[0] *= Yin/xyz[1];
+	xyz[2] *= Yin/xyz[1];
+	xyz[1] = Yin;
+}
+
 /* - - - - - - - - - - - - - - - - - - - - - - - - - */
 /* Spectral and illuminant chromaticity locus support */
 
@@ -5792,7 +6003,7 @@ static xslpoly illo_P_CIE_1931_2_uv       = { 2, icxOT_CIE_1931_2,     1, 0 };
 static xslpoly illo_P_CIE_1964_10_xy      = { 2, icxOT_CIE_1964_10,    0, 0 };
 static xslpoly illo_P_CIE_1964_10_uv      = { 2, icxOT_CIE_1964_10,    1, 0 };
 
-/* Return a pointer to the chromaticity locus poligon */
+/* Return a pointer to the (static) chromaticity locus poligon */
 /* return NULL on failure. */
 xslpoly *chrom_locus_poligon(
 icxLocusType loty,			/* Locus type, 1 = spectral, 2 = Daylight, 3 = Plankian */
@@ -6060,8 +6271,8 @@ double *in				/* Input XYZ values */
 /* Given an XYZ value, return sRGB values. */
 /* This is a little slow if wp used */
 void icx_XYZ2sRGB(
-double *out,			/* Return approximate sRGB values */
-double *wp,				/* Input XYZ white point (may be NULL) */
+double *out,			/* Return sRGB values */
+double *wp,				/* Input XYZ white point (D65 used if NULL) */
 double *in				/* Input XYZ values */
 ) {
 	int i, j;
@@ -6109,6 +6320,50 @@ double *in				/* Input XYZ values */
 	}
 }
 
+/* Given an RGB value, return XYZ values. */
+/* This is a little slow */
+void icx_sRGB2XYZ(
+double *out,			/* Return XYZ values */
+double *wp,				/* Output XYZ white point (D65 used if NULL, othewise Bradford) */
+double *in				/* Input sRGB values */
+) {
+	int i, j;
+	double tmp[3];
+	double d65[3] = { 0.950543, 1.0, 1.089303 };	/* D65 */
+	double imat[3][3] = {				/* sRGB absolute XYZ->RGB */
+		{ 0.4124, 0.3576, 0.1805 },
+		{ 0.2126, 0.7152, 0.0722 },
+		{ 0.0193, 0.1192, 0.9505 }
+	};
+
+	/* Undo gamma */
+	for (j = 0; j < 3; j++) {
+		if (in[j] < 0.04045)
+			tmp[j] = in[j]/12.92;
+		else
+			tmp[j] = pow((in[j] + 0.055)/1.055, 2.4);
+	}
+
+	/* Convert to XYZ cromaticities */
+	for (i = 0; i < 3; i++) {
+		out[i] = 0.0;
+		for (j = 0; j < 3; j++) {
+			out[i] += tmp[j] * imat[i][j];
+		}
+	}
+
+	/* Do a simple Bradford between D65 and wp */
+	if (wp != NULL) {
+		icmXYZNumber dst, src;
+		double vkmat[3][3];
+
+		icmAry2XYZ(src, d65);
+		icmAry2XYZ(dst, wp);
+		icmChromAdaptMatrix(ICM_CAM_BRADFORD, dst, src, vkmat);
+		icmMulBy3x3(out, vkmat, out);
+	}
+}
+
 /* Given an XYZ value, return approximate RGB value */
 /* Desaurate to white by the given amount */
 void icx_XYZ2RGB_ds(
diff --git a/xicc/xspect.h b/xicc/xspect.h
index e3adc1f..75c98f3 100644
--- a/xicc/xspect.h
+++ b/xicc/xspect.h
@@ -47,7 +47,7 @@ typedef struct {
 	int    spec_n;					/* Number of spectral bands, 0 if not valid */
 	double spec_wl_short;			/* First reading wavelength in nm (shortest) */
 	double spec_wl_long;			/* Last reading wavelength in nm (longest) */
-	double norm;					/* Normalising scale value */
+	double norm;					/* Normalising scale value, ie. 1, 100 etc. */
 	double spec[XSPECT_MAX_BANDS];	/* Spectral value, shortest to longest */
 } xspect;
 
@@ -113,6 +113,9 @@ void xspect_denorm(xspect *sp);
 #ifndef SALONEINSTLIB
 /* Convert from one xspect type to another */
 void xspect2xspect(xspect *dst, xspect *targ, xspect *src);
+
+/* Plot up to 3 spectra */
+void xspect_plot(xspect *sp1, xspect *sp2, xspect *sp3);
 #endif /* !SALONEINSTLIB*/
 
 /* ------------------------------------------------------------------------------ */
@@ -129,16 +132,18 @@ typedef enum {
 	icxIT_C          = 4,	/* Standard Illuminant C */
     icxIT_D50		 = 5,	/* Daylight 5000K */
     icxIT_D50M2		 = 6,	/* Daylight 5000K, UV filtered (M2) */
-    icxIT_D65		 = 7,	/* Daylight 6500K */
-    icxIT_E		     = 8,	/* Equal Energy */
+    icxIT_D55		 = 7,	/* Daylight 5500K (use specified temperature) */
+    icxIT_D65		 = 8,	/* Daylight 6500K */
+    icxIT_D75		 = 9,	/* Daylight 7500K (uses specified temperature) */
+    icxIT_E		     = 10,	/* Equal Energy = flat = 1.0 */
 #ifndef SALONEINSTLIB
-    icxIT_F5		 = 9,	/* Fluorescent, Standard, 6350K, CRI 72 */
-    icxIT_F8		 = 10,	/* Fluorescent, Broad Band 5000K, CRI 95 */
-    icxIT_F10		 = 11,	/* Fluorescent Narrow Band 5000K, CRI 81 */
-	icxIT_Spectrocam = 12,	/* Spectrocam Xenon Lamp */
-    icxIT_Dtemp		 = 13,	/* Daylight at specified temperature */
+    icxIT_F5		 = 11,	/* Fluorescent, Standard, 6350K, CRI 72 */
+    icxIT_F8		 = 12,	/* Fluorescent, Broad Band 5000K, CRI 95 */
+    icxIT_F10		 = 13,	/* Fluorescent Narrow Band 5000K, CRI 81 */
+	icxIT_Spectrocam = 14,	/* Spectrocam Xenon Lamp */
+    icxIT_Dtemp		 = 15,	/* Daylight at specified temperature */
 #endif /* !SALONEINSTLIB*/
-    icxIT_Ptemp		 = 14	/* Planckian at specified temperature */
+    icxIT_Ptemp		 = 16	/* Planckian at specified temperature */
 } icxIllumeType;
 
 /* Fill in an xpsect with a standard illuminant spectrum */
@@ -200,6 +205,7 @@ struct _xsp2cie {
 	xspect emits;	/* Estimated FWA emmission spectrum */
 	xspect media;	/* Estimated base media (ie. minus FWA) */
 	xspect tillum;	/* Y = 1 Normalised target/simulated instrument illuminant spectrum */
+					/* Use oillum if tillum spec_n = 0 */
 	xspect oillum;	/* Y = 1 Normalised observer illuminant spectrum */
 	double Sm;		/* FWA Stimulation level for emits contribution */
 	double FWAc;	/* FWA content (informational) */
@@ -230,14 +236,13 @@ struct _xsp2cie {
 	                 xspect *in				/* Spectrum to be converted, normalised by norm */
 	                );
 
-#ifndef SALONEINSTLIB
-	/* Set Media White. This enables extracting and applying the */
-	/* colorant reflectance value from/to the meadia. */
-	/* return NZ if error */
-	int (*set_mw) (struct _xsp2cie *p,	/* this */
-	                xspect *white		/* Spectrum of plain media */
-	                );
+	/* Get the XYZ of the illuminant being used to compute the CIE XYZ */
+	/* value. */
+	void (*get_cie_il)(struct _xsp2cie *p,	/* this */
+	                   double *xyz			/* Return the XYZ */
+	                   );
 
+#ifndef SALONEINSTLIB
 	/* Set Fluorescent Whitening Agent compensation */
 	/* return NZ if error */
 	int (*set_fwa) (struct _xsp2cie *p,	/* this */
@@ -247,13 +252,14 @@ struct _xsp2cie {
 	                xspect *white		/* Spectrum of plain media */
 	                );
 
-	/* Set FWA given updated conversion illuminant. */
+	/* Set FWA given updated conversion illuminants. */
 	/* (We assume that xsp2cie_set_fwa has been called first) */
 	/* return NZ if error */
 	int (*update_fwa_custillum) (struct _xsp2cie *p,
 					xspect *tillum,		/* Spectrum of target/simulated instrument illuminant, */
 										/* NULL to use set_fwa() value. */
 	                xspect *custIllum	/* Spectrum of observer illuminant */
+										/* NULL to use new_xsp2cie() value. */
 	                );	
 
 	/* Get Fluorescent Whitening Agent compensation information */
@@ -262,6 +268,14 @@ struct _xsp2cie {
 					double *FWAc		/* FWA content as a ratio. */
 	                );
 
+
+	/* Set Media White. This enables extracting and applying the */
+	/* colorant reflectance value from/to the meadia. */
+	/* return NZ if error */
+	int (*set_mw) (struct _xsp2cie *p,	/* this */
+	                xspect *white		/* Spectrum of plain media */
+	                );
+
 	/* Extract the colorant reflectance value from the media. Takes FWA */
 	/* into account if set. Media white or FWA must be set. */
 	/* return NZ if error */
@@ -309,6 +323,20 @@ double xyz[3],			/* Input XYZ value */
 int viscct				/* nz to use visual CIEDE2000, 0 to use CCT CIE 1960 UCS. */
 );
 
+/* Given a choice of temperature dependent illuminant (icxIT_Dtemp or icxIT_Ptemp), */
+/* a color temperature and a Y value, return the corresponding XYZ */
+/* An observer type can be chosen for interpretting the spectrum of the input and */
+/* the illuminant. */
+/* Return xyz[0] = -1.0 on erorr */
+void icx_ill_ct2XYZ(
+double xyz[3],			/* Return the XYZ value */
+icxIllumeType ilType,	/* Type of illuminant, icxIT_Dtemp or icxIT_Ptemp */
+icxObserverType obType,	/* Observer, CIE_1931_2 or CIE_1964_10 */
+int viscct,				/* nz to use visual CIEDE2000, 0 to use CCT CIE 1960 UCS. */
+double tin,				/* Input temperature */
+double Yin				/* Input Y value */
+);
+
 /* --------------------------- */
 /* Spectrum locus              */
 
@@ -334,7 +362,7 @@ typedef enum {
     icxLT_plankian	 = 3
 } icxLocusType;
 
-/* Return a pointer to the chromaticity locus object */
+/* Return a pointer to the (static) chromaticity locus object */
 /* return NULL on failure. */
 xslpoly *chrom_locus_poligon(icxLocusType locus_type, icxObserverType obType, int cspace);
 
@@ -371,10 +399,18 @@ double *in				/* Input XYZ values */
 /* return sRGB values */
 void icx_XYZ2sRGB(
 double *out,			/* Return sRGB value */
-double *wp,				/* Input XYZ white point (may be NULL) */
+double *wp,				/* Input XYZ white point (D65 used if NULL) */
 double *in				/* Input XYZ values */
 );
 
+/* Given an RGB value, return XYZ values. */
+/* This is a little slow */
+void icx_sRGB2XYZ(
+double *out,			/* Return XYZ values */
+double *wp,				/* Output XYZ white point (D65 used if NULL) */
+double *in				/* Input sRGB values */
+);
+
 /* Given an XYZ value, return approximate RGB value */
 /* Desaurate to white by the given amount */
 void icx_XYZ2RGB_ds(