New release 1.7.0

16 files changed, 1314 insertions, 845 deletions

diff --git a/profile/Jamfile b/profile/Jamfile
index 5c286b9..0c857a8 100644
--- a/profile/Jamfile
+++ b/profile/Jamfile
@@ -6,7 +6,7 @@ PREF_LINKFLAGS	+= $(LINKDEBUGFLAG) ;
 
 #Products
 Libraries = libprof ;
-Executables = cb2ti3 kodak2ti3 txt2ti3 splitti3 mppcheck mppprof
+Executables = cb2ti3 kodak2ti3 txt2ti3 ls2ti3 splitti3 mppcheck mppprof
               profcheck invprofcheck colverify colprof printcal applycal ;
 Headers = prof.h ;
 Samples = example.sp example121.sp 3dap5k.sp GTIPlus.sp Office.sp Trulux.sp TruluxPlus.sp
@@ -28,7 +28,7 @@ Library libprof : profin.c profout.c ;
 
 
 LINKLIBS = ../rspl/librspl ../icc/libicc ../cgats/libcgats ../numlib/libnum ../plot/libplot
-           ../plot/libvrml ;
+           ../plot/libvrml ../numlib/libui ;
 
 # Simple profile generator
 Main simpprof : simpprof.c ;
@@ -40,11 +40,14 @@ Main kodak2ti3 : kodak2ti3.c ;
 Main cb2ti3 : cb2ti3.c ;
 
 # the gcc linker is retarded, and can't link to things it's gone past, hence 2 x libxicc...
-LINKLIBS = ../xicc/libxicc ../spectro/libinsttypes ../xicc/libxicc ../gamut/libgamut $(LINKLIBS) ;
+LINKLIBS = ../xicc/libxicc ../spectro/libinsttypes ../spectro/libdisptechs ../xicc/libxicc ../gamut/libgamut $(LINKLIBS) ;
 
 #Gretag/Logo raw CMYK profile data to Argyll CGATS format converter
 Main txt2ti3 : txt2ti3.c ;
 
+#LightSpace to Argyll CGATS format
+Main ls2ti3 : ls2ti3.c : : : ../xml : : ../xml/libmxml ;
+
 #Split a .ti3 into two pieces randomly
 Main splitti3 : splitti3.c ;
 
@@ -85,7 +88,7 @@ if $(HOME) = "d:\\usr\\graeme" && $(PWD) = "/src/argyll/profile" {
 }
 
 # Development code
-if [ GLOB . : retargti3.c ]  {
+if [ GLOB [ NormPaths . ] : retargti3.c ]  {
 	Main retargti3 : retargti3.c ;
 }
 
diff --git a/profile/afiles b/profile/afiles
index 8b3b980..2d91df7 100644
--- a/profile/afiles
+++ b/profile/afiles
@@ -15,6 +15,7 @@ mppprof.c
 mppcheck.c
 simpprof.c
 kodak2ti3.c
+ls2ti3.c
 cb2ti3.c
 txt2ti3.c
 splitti3.c
diff --git a/profile/applycal.c b/profile/applycal.c
index ff9a76d..17187db 100644
--- a/profile/applycal.c
+++ b/profile/applycal.c
@@ -33,6 +33,7 @@
 #include "numlib.h"
 #include "rspl.h"
 #include "xicc.h"
+#include "ui.h"
 
 #undef DEBUG
 
diff --git a/profile/colprof.c b/profile/colprof.c
index b80f13c..1e73cdf 100644
--- a/profile/colprof.c
+++ b/profile/colprof.c
@@ -55,6 +55,7 @@
 #include "cgats.h"
 #include "xicc.h"
 #include "prof.h"
+#include "ui.h"
 
 #define DEFAVGDEV 0.5		/* Default average deviation percentage */
 							/* This equates to a uniform added error of +/- 1% */
@@ -66,7 +67,7 @@
   Flags used:
 
          ABCDEFGHIJKLMNOPQRSTUVWXYZ
-  upper  . ..    . ... .. ....    .
+  upper  ....    . ... .. ....    .
   lower  .... .. . .. .........    
 
 */
@@ -100,7 +101,7 @@ void usage(char *diag, ...) {
 	fprintf(stderr," -np             Don't create input (Device) grid position curves\n");
 	fprintf(stderr," -no             Don't create output (PCS) shaper curves\n");
 	fprintf(stderr," -nc             Don't put the input .ti3 data in the profile\n");
-	fprintf(stderr," -k zhxr         Black value target: z = zero K,\n");
+	fprintf(stderr," -k zhxr         Black Ink generation target: z = zero K,\n");
 	fprintf(stderr,"                 h = 0.5 K, x = max K, r = ramp K (def.)\n");
 	fprintf(stderr," -k p stle stpo enpo enle shape\n");
 	fprintf(stderr,"                 stle: K level at White 0.0 - 1.0\n");
@@ -122,6 +123,7 @@ void usage(char *diag, ...) {
 	fprintf(stderr," -uc             If input profile, clip cLUT values above WP\n");
 	fprintf(stderr," -U scale        If input profile, scale media white point by scale\n");
 	fprintf(stderr," -R              Restrict white <= 1.0, black and primaries to be +ve\n");
+//	fprintf(stderr," -B X,Y,Z        Display Black Point override hack\n");
 	fprintf(stderr," -V demphasis    Degree of dark region cLUT grid emphasis 1.0-4.0 (default %.2f = none)\n",DEMPH_DEFAULT);
 	fprintf(stderr," -f [illum]      Use Fluorescent Whitening Agent compensation [opt. simulated inst. illum.:\n");
 	fprintf(stderr,"                  M0, M1, M2, A, C, D50 (def.), D50M2, D65, F5, F8, F10 or file.sp]\n");
@@ -185,6 +187,7 @@ int main(int argc, char *argv[]) {
 	int autowpsc = 0;			/* Auto scale the WP to prevent clipping above WP patch */
 	int clipovwp = 0;			/* Clip cLUT values above WP */
 	int clipprims = 0;			/* Clip white, black and primaries */
+//	double bpo[3] = { -1,-1,-1 };	/* Black point override hack XYZ value */
 	double demph = 0.0;			/* Emphasise dark region grid resolution in cLUT */
 	double iwpscale = -1.0;		/* Input white point scale factor */
 	int doinextrap = 1;			/* Sythesize extra sample points for input device cLUT */
@@ -200,6 +203,7 @@ int main(int argc, char *argv[]) {
 	int spec = 0;				/* Use spectral data flag */
 	icxIllumeType tillum = icxIT_none;	/* Target/simulated instrument illuminant */ 
 	xspect cust_tillum;			/* Custom target/simulated illumination spectrum */
+								/* xspect will use illum/cust_illum if tillum == none */
 	icxIllumeType illum = icxIT_D50;	/* Spectral defaults */
 	xspect cust_illum;			/* Custom illumination spectrum */
 	icxObserverType observ = icxOT_CIE_1931_2;	/* The classic observer */
@@ -428,11 +432,6 @@ int main(int argc, char *argv[]) {
 					oquality = 0;
 			}
 
-			else if (argv[fa][1] == 'B') {
-				oquality = -2;
-				doinb2a = 0;
-			}
-
 			/* Disable input or output luts */
 			else if (argv[fa][1] == 'n') {
 				fa = nfa;
@@ -478,6 +477,18 @@ int main(int argc, char *argv[]) {
 				clipprims = 1;
 			}
 
+#ifdef NEVER	/* Prototype - not used */
+			/* Black Point override hack */
+			else if (argv[fa][1] == 'B') {
+				if (na == NULL) usage("Expect X,Y,Z value after -B");
+				fa = nfa;
+				if (sscanf(na, " %lf , %lf , %lf ",&bpo[0],&bpo[1],&bpo[2]) != 3)
+					usage("Couldn't parse hack black point (-B) value '%s'",na);
+				if (bpo[0] < 0.0 || bpo[1] < 0.0 || bpo[1] < 0.0)
+					usage("Bad hack black point (-B) value '%s'",na);
+			}
+#endif
+
 			/* Degree of dark region emphasis */
 			else if (argv[fa][1] == 'V') {
 				if (na == NULL) usage(0,"Expected argument to dark emphasis flag -V");
@@ -824,7 +835,7 @@ int main(int argc, char *argv[]) {
 					if (sscanf(na+1,":%lf:%lf:%lf",&x,&y,&z) == 3) {
 						vc->Wxyz[0] = x; vc->Wxyz[1] = y; vc->Wxyz[2] = z;
 					} else if (sscanf(na+1,":%lf:%lf",&x,&y) == 2) {
-						vc->Wxyz[0] = x; vc->Wxyz[1] = y;
+						vc->Wxyz[0] = x; vc->Wxyz[1] = y; vc->Wxyz[2] = -1;
 					} else
 						usage("Viewing condition (-%cw) unrecognised white point '%s'",argv[fa][1],na+1);
 				} else if (na[0] == 'a' || na[0] == 'A') {
@@ -842,13 +853,13 @@ int main(int argc, char *argv[]) {
 				} else if (na[0] == 'f' || na[0] == 'F') {
 					if (na[1] != ':')
 						usage("Viewing conditions (-%cf) missing ':'",argv[fa][1]);
-					vc->Yf = atof(na+2);
+					vc->Yf = atof(na+2)/100.0;
 				} else if (na[0] == 'g' || na[0] == 'G') {
 					double x, y, z;
 					if (sscanf(na+1,":%lf:%lf:%lf",&x,&y,&z) == 3) {
 						vc->Gxyz[0] = x; vc->Gxyz[1] = y; vc->Gxyz[2] = z;
 					} else if (sscanf(na+1,":%lf:%lf",&x,&y) == 2) {
-						vc->Gxyz[0] = x; vc->Gxyz[1] = y;
+						vc->Gxyz[0] = x; vc->Gxyz[1] = y; vc->Gxyz[2] = -1;
 					} else if (sscanf(na+1,":%lf",&x) == 1) {
 						vc->Yg = x/100.0;
 					} else
@@ -966,6 +977,9 @@ int main(int argc, char *argv[]) {
 	if (strcmp(icg->t[0].kdata[ti],"OUTPUT") == 0) {
 		icxInk ink;							/* Ink parameters */
 
+//		if (bpo[1] >= 0.0)
+//			error("-B option not valid for output profile");
+
 		if ((ti = icg->find_kword(icg, 0, "TOTAL_INK_LIMIT")) >= 0) {
 			int imax;
 			imax = atoi(icg->t[0].kdata[ti]);
@@ -1075,7 +1089,9 @@ int main(int argc, char *argv[]) {
 
 		make_output_icc(ptype, 0, iccver, verb, iquality, oquality,
 		                noisluts, noipluts, nooluts, nocied, noptop, nostos,
-		                gamdiag, verify, clipprims, iwpscale, &ink, inname, outname, icg, 
+		                gamdiag, verify, clipprims, iwpscale,
+//		                NULL,		/* bpo */
+		                &ink, inname, outname, icg, 
 		                spec, tillum, &cust_tillum, illum, &cust_illum, observ, fwacomp,
 		                smooth, avgdev, 1.0,
 		                ipname[0] != '\000' ? ipname : NULL,
@@ -1086,6 +1102,9 @@ int main(int argc, char *argv[]) {
 
 	} else if (strcmp(icg->t[0].kdata[ti],"INPUT") == 0) {
 
+//		if (bpo[1] >= 0.0)
+//			error("-B option not valid for input profile");
+
 		if (ptype == prof_default)
 			ptype = prof_clutLab;		/* For best possible quality */
 
@@ -1113,7 +1132,9 @@ int main(int argc, char *argv[]) {
 		/* If a source gamut is provided for a Display, then a V2.4.0 profile will be created */
 		make_output_icc(ptype, mtxtoo, iccver, verb, iquality, oquality,
 		                noisluts, noipluts, nooluts, nocied, noptop, nostos,
-		                gamdiag, verify, clipprims, iwpscale, NULL, inname, outname, icg,
+		                gamdiag, verify, clipprims, iwpscale,
+//		                bpo[1] >= 0.0 ? bpo : NULL,
+		                NULL, inname, outname, icg,
 		                spec, icxIT_none, NULL, illum, &cust_illum, observ, 0,
 		                smooth, avgdev, demph,
 		                ipname[0] != '\000' ? ipname : NULL,
diff --git a/profile/colverify.c b/profile/colverify.c
index ba09366..e84fa6a 100644
--- a/profile/colverify.c
+++ b/profile/colverify.c
@@ -22,7 +22,7 @@
  * TTBD:
  */
 
-#undef DEBUG
+#define DEBUG
 
 #define verbo stdout
 
@@ -37,9 +37,20 @@
 #include "vrml.h"
 #include "cgats.h"
 #include "xicc.h"
-#include "ccmx.h"
 #include "insttypes.h"
+#include "disptechs.h"
+#include "ccmx.h"
 #include "sort.h"
+#include "plot.h"
+#include "ui.h"
+
+#ifdef DEBUG
+#undef DBG
+#define DBG(xxx) printf xxx ;
+#else
+#undef DBG
+#define DBG(xxx) 
+#endif
 
 void
 usage(void) {
@@ -50,13 +61,17 @@ usage(void) {
 	fprintf(stderr," -n              Normalise each files reading to its white Y\n");
 	fprintf(stderr," -N              Normalise each files reading to its white XYZ\n");
 	fprintf(stderr," -m              Normalise each files reading to its white X+Y+Z\n");
+	fprintf(stderr," -M              Normalise both files reading to mean white XYZ\n");
 	fprintf(stderr," -D              Use D50 100.0 as L*a*b* white reference\n");
 	fprintf(stderr," -c              Show CIE94 delta E values\n");
 	fprintf(stderr," -k              Show CIEDE2000 delta E values\n");
+	fprintf(stderr," -h              Plot a histogram of delta E's\n");
 	fprintf(stderr," -s              Sort patch values by error\n");
-	fprintf(stderr," -w              create VRML vector visualisation (measured.wrl)\n");
-	fprintf(stderr," -W              create VRML marker visualisation (measured.wrl)\n");
-	fprintf(stderr," -x              Use VRML axes\n");
+	fprintf(stderr," -w              create PCS %s vector visualisation (measured%s)\n",vrml_format(), vrml_ext());
+	fprintf(stderr," -W              create PCS %s marker visualisation (measured%s)\n",vrml_format(),vrml_ext());
+	fprintf(stderr," -d              create Device RGB %s marker visualisation (measured%s)\n",vrml_format(),vrml_ext());
+//	fprintf(stderr," -d y            create Device YCbCr %s marker visualisation (measured%s)\n",vrml_format(),vrml_ext());
+	fprintf(stderr," -x              Use %s axes\n",vrml_format());
 	fprintf(stderr," -f [illum]      Use Fluorescent Whitening Agent compensation [opt. simulated inst. illum.:\n");
 	fprintf(stderr,"                  M0, M1, M2, A, C, D50 (def.), D50M2, D65, F5, F8, F10 or file.sp]\n");
 	fprintf(stderr," -i illum        Choose illuminant for computation of CIE XYZ from spectral data & FWA:\n");
@@ -74,6 +89,8 @@ usage(void) {
 typedef struct {
 	char sid[50];		/* sample id */
 	char loc[100];		/* sample location (empty if none) */
+	double rgb[3];		/* RGB value if RGB device space present, or YCbCr if dovrml==4 */
+	double ycc[3];		/* YCbCr if RGB and dovrml==4 */
 	int og;				/* Out of gamut flag */
 	double xyz[3];		/* XYZ value */
 	double v[3];		/* Lab value */
@@ -82,17 +99,25 @@ typedef struct {
 	double ide[3];		/* Lab Component DE */
 } pval;
 
+/* Histogram bin type */
+typedef struct {
+	int count;			/* Raw count */
+	double val;			/* Normalized value */
+	double min, max;	/* Bin range */
+} hbin;
+
 int main(int argc, char *argv[])
 {
 	int fa,nfa,mfa;			/* current argument we're looking at */
 	int verb = 0;       	/* Verbose level */
 	int norm = 0;			/* 1 = norm to White Y, 2 = norm to White XYZ */
-							/* 3 = norm to White X+Y+Z */
+							/* 3 = norm to White X+Y+Z, 4 = norm to average XYZ */
 	int usestdd50 = 0;		/* Use standard D50 instead of avg white as reference */
 	int cie94 = 0;
 	int cie2k = 0;
-	int dovrml = 0;
+	int dovrml = 0;			/* 1 = PCS vector, 2 = PCS marker, 3 = RGB, 4 - YCbCr */
 	int doaxes = 0;
+	int dohisto = 0;			/* Plot histogram of delta E's */
 	int dosort = 0;
 	char ccmxname[MAXNAMEL+1] = "\000";  /* Colorimeter Correction Matrix name */
 	ccmx *cmx = NULL;					/* Colorimeter Correction Matrix */
@@ -106,6 +131,7 @@ int main(int argc, char *argv[])
 		char name[MAXNAMEL+1];	/* Patch filename  */
 		int isdisp;				/* nz if display */
 		int isdnormed;      	/* Has display data been normalised to 100 ? */
+		int isrgb;				/* Is RGB device space ? */
 		int npat;				/* Number of patches */
 		int nig;				/* Number of patches in gamut */
 		double w[3];			/* XYZ of "white" */
@@ -125,16 +151,11 @@ int main(int argc, char *argv[])
 
 	icmXYZNumber labw = icmD50;	/* The Lab white reference */
 
-	char out_name[MAXNAMEL+4+1]; /* VRML name */
+	char out_name[MAXNAMEL+4+1]; /* VRML/X3D name */
 	vrml *wrl = NULL;
 
 	int i, j, n;
 
-#if defined(__IBMC__)
-	_control87(EM_UNDERFLOW, EM_UNDERFLOW);
-	_control87(EM_OVERFLOW, EM_OVERFLOW);
-#endif
-
 	if (argc <= 1)
 		usage();
 
@@ -182,15 +203,31 @@ int main(int argc, char *argv[])
 				norm = 3;
 			}
 
+			else if (argv[fa][1] == 'M') {
+				norm = 4;
+			}
+
 			else if (argv[fa][1] == 'D')
 				usestdd50 = 1;
 
-			/* VRML */
+			/* VRML/X3D */
 			else if (argv[fa][1] == 'w')
 				dovrml = 1;
+
 			else if (argv[fa][1] == 'W')
 				dovrml = 2;
 
+			else if (argv[fa][1] == 'd') {
+				dovrml = 3;
+				if (na != NULL) {	/* Argument is present - RGB or YCbCr. */
+					fa = nfa;
+					if (strcmp(na, "y") == 0)
+						dovrml = 4;
+					else
+						usage();
+				}					
+			}
+
 			/* Axes */
 			else if (argv[fa][1] == 'x')
 				doaxes = 1;
@@ -206,6 +243,10 @@ int main(int argc, char *argv[])
 				cie2k = 1;
 			}
 
+			/* Plot histogram */
+			else if (argv[fa][1] == 'h')
+				dohisto = 1;
+
 			/* Sort */
 			else if (argv[fa][1] == 's')
 				dosort = 1;
@@ -337,13 +378,13 @@ int main(int argc, char *argv[])
 	if (fa >= argc || argv[fa][0] == '-') usage();
 	strncpy(cg[1].name,argv[fa],MAXNAMEL); cg[1].name[MAXNAMEL] = '\000';
 
-	/* Create VRML name */
+	/* Create VRML/X3D base name */
 	{
 		char *xl;
 		strcpy(out_name, cg[1].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 */
 	}
 
 	if (fwacomp && spec == 0)
@@ -396,11 +437,14 @@ int main(int argc, char *argv[])
 		int sidx;					/* Sample ID index */
 		int sldx = -1;				/* Sample location index, < 0 if invalid */
 		int xix, yix, zix;
+		int rgbix[3];				/* RGB field indexes (if rgb ) */
 
 		/* Open CIE target values */
 		cgf = new_cgats();			/* Create a CGATS structure */
 		cgf->add_other(cgf, ""); 	/* Allow any signature file */
 	
+		DBG(("Opening file '%s'\n",cg[n].name))
+
 		if (cgf->read_name(cgf, cg[n].name))
 			error("CGATS file '%s' read error : %s",cg[n].name,cgf->err);
 	
@@ -438,34 +482,50 @@ int main(int argc, char *argv[])
 			cg[n].isdnormed = 0;
 			cg[n].w[0] = cg[n].w[1] = cg[n].w[2] = 0.0;
 
-			if ((ti = cgf->find_kword(cgf, 0, "DEVICE_CLASS")) < 0)
-				error ("Input file '%s' doesn't contain keyword DEVICE_CLASS",cg[n].name);
-	
-			if (strcmp(cgf->t[0].kdata[ti],"DISPLAY") == 0) {
-				cg[n].isdisp = 1;
-				cg[n].isdnormed = 1;	/* Assume display type is normalised to 100 */
-				illum = icxIT_none;		/* Displays are assumed to be self luminous */
-				/* ?? What if two files are different ?? */
-			}
-
-			if (cg[n].isdisp) {
+			if ((ti = cgf->find_kword(cgf, 0, "DEVICE_CLASS")) < 0) {
+				warning("Input file '%s' doesn't contain keyword DEVICE_CLASS",cg[n].name);
 
-				if ((ti = cgf->find_kword(cgf, 0, "LUMINANCE_XYZ_CDM2")) >= 0) {
-					if (sscanf(cgf->t[0].kdata[ti], " %lf %lf %lf ",&cg[n].w[0], &cg[n].w[1], &cg[n].w[2]) != 3)
-						cg[n].w[0] = cg[n].w[1] = cg[n].w[2] = 0.0;
+			} else {
+				if (strcmp(cgf->t[0].kdata[ti],"DISPLAY") == 0) {
+					cg[n].isdisp = 1;
+					cg[n].isdnormed = 1;	/* Assume display type is normalised to 100 */
+					illum = icxIT_none;		/* Displays are assumed to be self luminous */
+					/* ?? What if two files are different ?? */
 				}
-
-				/* See if there is an explicit tag indicating data has been normalised to Y = 100 */
-				if ((ti = cgf->find_kword(cgf, 0, "NORMALIZED_TO_Y_100")) >= 0) {
-					if (strcmp(cgf->t[0].kdata[ti],"NO") == 0) {
-						cg[n].isdnormed = 0;
-					} else {
-						cg[n].isdnormed = 1;
+	
+				if (cg[n].isdisp) {
+	
+					if ((ti = cgf->find_kword(cgf, 0, "LUMINANCE_XYZ_CDM2")) >= 0) {
+						if (sscanf(cgf->t[0].kdata[ti], " %lf %lf %lf ",&cg[n].w[0], &cg[n].w[1], &cg[n].w[2]) != 3)
+							cg[n].w[0] = cg[n].w[1] = cg[n].w[2] = 0.0;
+					}
+	
+					/* See if there is an explicit tag indicating data has been normalised to Y = 100 */
+					if ((ti = cgf->find_kword(cgf, 0, "NORMALIZED_TO_Y_100")) >= 0) {
+						if (strcmp(cgf->t[0].kdata[ti],"NO") == 0) {
+							cg[n].isdnormed = 0;
+						} else {
+							cg[n].isdnormed = 1;
+						}
 					}
 				}
 			}
 		}
 
+		/* See if it has RGB device space (for -d option) */
+		if ((rgbix[0] = cgf->find_field(cgf, 0, "RGB_R")) >= 0
+		 && cgf->t[0].ftype[rgbix[0]] == r_t
+
+		 &&	(rgbix[1] = cgf->find_field(cgf, 0, "RGB_G")) >= 0
+		 && cgf->t[0].ftype[rgbix[1]] == r_t
+
+		 &&	(rgbix[2] = cgf->find_field(cgf, 0, "RGB_B")) >= 0
+		 && cgf->t[0].ftype[rgbix[2]] == r_t) {
+			cg[n].isrgb = 1;
+		} else {
+			cg[n].isrgb = 0;
+		}
+
 		/* Read all the target patches */
 		if (cg[n].npat <= 0)
 			error("No sets of data in file '%s'",cg[n].name);
@@ -535,14 +595,14 @@ int main(int argc, char *argv[])
 				cg[n].pat[i].xyz[1] = *((double *)cgf->t[0].fdata[i][yix]);
 				cg[n].pat[i].xyz[2] = *((double *)cgf->t[0].fdata[i][zix]);
 
-				if (isLab) {	/* Convert to XYZ */
+				if (isLab) {	/* Convert Lab to XYZ */
 					icmLab2XYZ(&icmD50, cg[n].pat[i].xyz, cg[n].pat[i].xyz);
 				}
 //printf("~1 file %d patch %d = XYZ %f %f %f\n", n,i,cg[n].pat[i].xyz[0],cg[n].pat[i].xyz[1],cg[n].pat[i].xyz[2]);
 
 				/* restore normalised display values to absolute */
 				if (cg[n].isdnormed) {
-					if (cg[n].w[1] > 0.0) {
+					if (cg[n].w[1] > 0.0) {		// Found absoluute display white tag
 						cg[n].pat[i].xyz[0] *= cg[n].w[1]/100.0;
 						cg[n].pat[i].xyz[1] *= cg[n].w[1]/100.0;
 						cg[n].pat[i].xyz[2] *= cg[n].w[1]/100.0;
@@ -560,6 +620,16 @@ int main(int argc, char *argv[])
 				if (n == 1 && cmx != NULL) {
 					cmx->xform(cmx, cg[n].pat[i].xyz, cg[n].pat[i].xyz);
 				}
+
+				if ((dovrml == 3 || dovrml == 4) && cg[n].isrgb) {
+					cg[n].pat[i].rgb[0] = 0.01 * *((double *)cgf->t[0].fdata[i][rgbix[0]]);
+					cg[n].pat[i].rgb[1] = 0.01 * *((double *)cgf->t[0].fdata[i][rgbix[1]]);
+					cg[n].pat[i].rgb[2] = 0.01 * *((double *)cgf->t[0].fdata[i][rgbix[2]]);
+
+					if (dovrml == 4) {
+						icmRec709_RGBd_2_YPbPr(cg[n].pat[i].ycc, cg[n].pat[i].rgb);
+					}
+				}
 			}
 
 		} else { 		/* Using spectral data */
@@ -677,7 +747,7 @@ int main(int argc, char *argv[])
 
 				/* restore normalised display values to absolute */
 				if (cg[n].isdnormed) {
-					if (cg[n].w[1] > 0.0) {
+					if (cg[n].w[1] > 0.0) {		// Found absoluute display white tag
 						cg[n].pat[i].xyz[0] *= cg[n].w[1];
 						cg[n].pat[i].xyz[1] *= cg[n].w[1];
 						cg[n].pat[i].xyz[2] *= cg[n].w[1];
@@ -697,10 +767,11 @@ int main(int argc, char *argv[])
 
 
 		/* Locate the patch with maximum Y, a possible white patch */
-		if (norm) {
+		/* in case we need it latter. */
+		{
 			int ii;
 
-			if (cg[n].w[1] == 0.0) {	/* No white patch */
+			if (cg[n].w[1] == 0.0) {	/* No display white patch tag */
 
 				/* Locate patch with biggest Y, assume it is white... */
 				for (i = 0; i < cg[n].npat; i++) {
@@ -709,15 +780,25 @@ int main(int argc, char *argv[])
 						ii = i;
 					}
 				}
-				if (verb) printf("File %d Chose patch %d as white, xyz %f %f %f\n",
+				if (verb) printf("File %d Chose patch %d as white, XYZ %f %f %f\n",
 				                       n, ii+1,cg[n].w[0],cg[n].w[1],cg[n].w[2]);
 			} else {
-				if (verb) printf("File %d White is from display luminance ref. xyz %f %f %f\n",
+				if (verb) printf("File %d White is from display luminance ref. XYZ %f %f %f\n",
 				                       n, cg[n].w[0],cg[n].w[1],cg[n].w[2]);
 			}
 			icmCpy3(cg[n].nw, cg[n].w);
 		}
+		cgf->del(cgf);		/* Clean up */
+	}	/* Next file */
 
+	if (norm == 4) {		/* Normalise to average of white XYZ of the two files */
+		icmBlend3(cg[0].w, cg[0].w, cg[1].w, 0.5);
+		icmCpy3(cg[1].w, cg[0].w);
+//		if (verb) printf("Average White XYZ %f %f %f\n",cg[0].w[0],cg[0].w[1],cg[0].w[2]);
+	}
+
+	/* For both files */
+	for (n = 0; n < 2; n++) {
 
 		/* Normalise this file to white = 1.0 or D50 */
 		if (norm) {
@@ -725,7 +806,9 @@ int main(int argc, char *argv[])
 
 			double chmat[3][3];				/* Chromatic adapation matrix */
 
-			if (norm == 2) {		/* Norm to white XYZ */ 
+			DBG(("Normalizng '%s' to white\n",cg[n].name))
+
+			if (norm == 2 || norm == 4) {		/* Norm to white XYZ */ 
 				icmXYZNumber s_wp;
 				icmAry2XYZ(s_wp, cg[n].w);
 				icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, s_wp, chmat);
@@ -736,7 +819,7 @@ int main(int argc, char *argv[])
 					cg[n].pat[i].xyz[0] *= 100.0 / cg[n].w[1];
 					cg[n].pat[i].xyz[1] *= 100.0 / cg[n].w[1];
 					cg[n].pat[i].xyz[2] *= 100.0 / cg[n].w[1];
-				} else if (norm == 2) { 
+				} else if (norm == 2 || norm == 4) { 
 					icmMulBy3x3(cg[n].pat[i].xyz, chmat, cg[n].pat[i].xyz);
 				} else {
 					cg[n].pat[i].xyz[0] *= 100.0 / (cg[n].w[0] + cg[n].w[1] + cg[n].w[2]);
@@ -750,7 +833,7 @@ int main(int argc, char *argv[])
 				cg[n].nw[0] *= 100.0 / cg[n].w[1];
 				cg[n].nw[1] *= 100.0 / cg[n].w[1];
 				cg[n].nw[2] *= 100.0 / cg[n].w[1];
-			} else if (norm == 2) { 
+			} else if (norm == 2 || norm == 4) { 
 				icmMulBy3x3(cg[n].nw, chmat, cg[n].w);
 			} else {
 				cg[n].nw[0] *= 100.0 / (cg[n].w[0] + cg[n].w[1] + cg[n].w[2]);
@@ -759,8 +842,8 @@ int main(int argc, char *argv[])
 			}
 //printf("~1 file %d norm white XYZ %f %f %f\n", n,cg[n].nw[0], cg[n].nw[1], cg[n].nw[2]);
 		}
-		cgf->del(cgf);		/* Clean up */
-	}
+	}	/* Next file */
+
 	if (cmx != NULL) 
 		cmx->del(cmx);
 	cmx = NULL;
@@ -791,12 +874,24 @@ int main(int argc, char *argv[])
 		icmXYZNumber s_wp;
 		int rv;
 
+		DBG(("Figuring out of gamut patches\n"))
+
 		/* Convert sample PCS to relative */
-		icmAry2XYZ(s_wp, cg[0].nw);
+//printf("   cg[0].w %f %f %f\n", cg[0].w[0], cg[0].w[1], cg[0].w[2]);
+		icmAry2XYZ(s_wp, cg[0].w);
+		s_wp.X /= s_wp.Y;		// Normalise the white to 1.0
+		s_wp.Y /= s_wp.Y;		// so that matrix doesn't change magnitude
+		s_wp.Z /= s_wp.Y;
+//printf("   s_wp %f %f %f\n", s_wp.X, s_wp.Y, s_wp.Z);
 		icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, s_wp, chmat);
+//printf("~1 matrix = \n");
+//printf("   %f %f %f\n", chmat[0][0], chmat[0][1], chmat[0][2]);
+//printf("   %f %f %f\n", chmat[1][0], chmat[1][1], chmat[1][2]);
+//printf("   %f %f %f\n", chmat[2][0], chmat[2][1], chmat[2][2]);
 
 		for (i = 0; i < cg[0].npat; i++) {
 			icmMulBy3x3(in, chmat, cg[0].pat[i].xyz);
+
 //printf("~1 %d: xyz %f %f %f, rel %f %f %f\n", i+1, cg[0].pat[i].xyz[0], cg[0].pat[i].xyz[1], cg[0].pat[i].xyz[2], in[0], in[1], in[2]);
 
 			if ((rv = luo->inv_lookup(luo, out, in)) > 0 || 1) {
@@ -816,7 +911,13 @@ int main(int argc, char *argv[])
 			}
 		}
 		if (verb)
-			fprintf(verbo,"No of test patches in gamut = %d/%d\n",cg[0].npat - cg[0].nig,cg[0].npat);
+			fprintf(verbo,"No of test patches in gamut = %d/%d\n",cg[0].nig,cg[0].npat);
+	}
+
+	if (cg[0].nig <= 0) {
+		if (verb)
+			fprintf(verbo,"No test patches in gamut - givig up\n");
+		return 0;
 	}
 
 	/* Adjust the Lab reference white to be the mean of the white of the two files */
@@ -829,6 +930,8 @@ int main(int argc, char *argv[])
 			printf("L*a*b* white reference = XYZ %f %f %f\n",labw.X,labw.Y,labw.Z);
 	}
 
+	/* labw defaults to D50 */
+
 	/* Convert XYZ to Lab */
 	for (n = 0; n < 2; n++) {
 		for (i = 0; i < cg[n].npat; i++) {
@@ -837,11 +940,9 @@ int main(int argc, char *argv[])
 	}
 
 	/* Compute the delta E's */
+	DBG(("Computing the delta E's\n"))
 	for (i = 0; i < cg[0].npat; i++) {
 
-		if (cg[0].pat[i].og)		/* Skip out of gamut patches */
-			continue;
-
 		cg[0].pat[i].ixde[0] = fabs(cg[0].pat[i].xyz[0] - cg[1].pat[match[i]].xyz[0]);
 		cg[0].pat[i].ixde[1] = fabs(cg[0].pat[i].xyz[1] - cg[1].pat[match[i]].xyz[1]);
 		cg[0].pat[i].ixde[2] = fabs(cg[0].pat[i].xyz[2] - cg[1].pat[match[i]].xyz[2]);
@@ -869,6 +970,111 @@ int main(int argc, char *argv[])
 #undef HEAP_COMPARE
 
 	/* - - - - - - - - - - */
+	/* Plot a dE histogram */
+	if (dohisto) {
+		double demax = -1e6, demin = 1e6;
+		int maxbins = 50;		/* Maximum bins */
+		int minbins = 20;		/* Target minimum bins (depends on distribution) */ 
+		int mincount = 10;		/* Minimum number of points in a bin */
+		double mbwidth;
+		int nbins = 0;
+		hbin *bins = NULL;
+		pval **stpat;          /* Pointers to sorted cg[0].pat[] */
+		double tval;
+		double *x, *y;
+		
+		DBG(("Plotting histogram\n"))
+
+		/* Figure out the range of dE's */
+		for (i = 0; i < cg[0].npat; i++) {
+			double de = cg[0].pat[i].de;
+
+			if (de > demax)
+				demax = de;
+			if (de < demin)
+				demin = de;
+		}
+
+		if (demax < 1e-6)
+			error("histogram: dE range is too small to plot");
+
+		/* Bin width that gives maxbins */
+		mbwidth = demax / maxbins;
+		
+		/* Reduce mincount if needed to get minbins */
+		if (cg[0].npat/minbins < mincount)
+			mincount = cg[0].npat/minbins;
+
+		if ((bins = (hbin *)calloc(maxbins, sizeof(hbin))) == NULL)
+			error("malloc of histogram bins failed");
+
+		if ((stpat = (pval **)malloc(sizeof(pval *) * cg[0].npat)) == NULL)
+			error("Malloc failed - stpat[]");
+
+		for (i = 0; i < cg[0].npat; i++)
+			stpat[i] = &cg[0].pat[i];
+
+	  	/* Sort the dE's */
+#define HEAP_COMPARE(A,B) (A->de < B->de)
+		HEAPSORT(pval *, stpat, cg[0].npat);
+#undef HEAP_COMPARE
+
+		/* Create bins and add points */
+		bins[0].min = 0.0;
+		for (nbins = i = 0; i < cg[0].npat && nbins < maxbins; i++) {
+			double de = stpat[i]->de;
+
+			/* Move on to next bin ? */
+			if (bins[nbins].count >= mincount
+			 && (de - bins[nbins].min) >= mbwidth) {
+				if (i > 0)
+					bins[nbins].max = 0.5 * (de + stpat[i-1]->de);
+				else
+					bins[nbins].max = de;
+				nbins++;
+				bins[nbins].min = bins[nbins-1].max; 
+			} 
+			bins[nbins].count++;
+			bins[nbins].max = de;
+		}
+		if (bins[nbins].count != 0)
+			nbins++;
+
+		/* Compute value */
+		tval = 0.0;
+		for (i = 0; i < nbins; i++) {
+			bins[i].val = bins[i].count/(bins[i].max - bins[i].min);
+			tval += bins[i].val;
+		}
+
+		tval /= 100.0;		/* Make it % */
+		for (i = 0; i < nbins; i++) {
+			bins[i].val /= tval;
+			if (verb) fprintf(verbo,"Bin %d, %f - %f, % 2.4f%%, count %d\n",
+			             i,bins[i].min,bins[i].max,bins[i].val,bins[i].count);
+		}
+
+		/* Plot it */
+		if ((x = (double *)calloc(nbins+1, sizeof(double))) == NULL)
+			error("malloc of histogram x array");
+		if ((y = (double *)calloc(nbins+1, sizeof(double))) == NULL)
+			error("malloc of histogram y array");
+
+		for (i = 0; i < nbins; i++) {
+			x[i] = 0.5 * (bins[i].min + bins[i].max);
+			y[i] = bins[i].val;
+		}
+		x[i] = demax;
+		y[i] = 0.0;
+		do_plot(x, y, NULL, NULL, nbins+1);
+
+		free(y);
+		free(x);
+		free(bins);
+		free(stpat);
+	}
+
+	/* - - - - - - - - - - */
 	/* Figure out the report */
 	{
 		double merr = 0.0, aerr = 0.0;
@@ -879,13 +1085,60 @@ int main(int argc, char *argv[])
 		double rad;
 		double aierr[3] = { 0.0, 0.0, 0.0 };
 		double aixerr[3] = { 0.0, 0.0, 0.0 };
+		double red[3] = { 1.0, 0.2, 0.2 };
+		double green[3] = { 0.2, 1.0, 0.2 };
+		double min[3], max[3];
+		double col[3];
 
 		if (dovrml) {
-			wrl = new_vrml(out_name, doaxes, 0);
+			double vol;
+			int k;
+
+			wrl = new_vrml(out_name, doaxes, (dovrml == 3 || dovrml == 4) ? vrml_rgb : vrml_lab);
 			wrl->start_line_set(wrl, 0);
 
-			/* Fudge sphere diameter */
-			rad = 10.0/pow(cg[0].npat, 1.0/3.0);
+			for (j = 0; j < 3; j++) {
+				min[j] = 1e6;
+				max[j] = -1e6;
+			}
+
+			/* Get bounding box */
+			for (i = 0; i < cg[0].npat; i++) {
+				for (k = 0; k < 2; k++) {
+					for (j = 0; j < 3; j++) {
+						if (dovrml == 3 || dovrml == 4) {		/* RGB or YCC device plot */
+							if (cg[k].pat[i].rgb[j] > max[j])
+								max[j] = cg[k].pat[i].rgb[j];
+							if (cg[k].pat[i].rgb[j] < min[j])
+								min[j] = cg[k].pat[i].rgb[j];
+						} else {
+							if (cg[k].pat[i].v[j] > max[j])
+								max[j] = cg[k].pat[i].v[j];
+							if (cg[k].pat[i].v[j] < min[j])
+								min[j] = cg[k].pat[i].v[j];
+						}
+					}
+				}
+			}
+
+			for (vol = 1.0, j = 0; j < 3; j++) {
+//printf("~1 size[%d] = %f\n",j, max[j] - min[j]);
+				vol *= (max[j] - min[j]);
+			}
+			vol = sqrt(vol);
+//printf("~1 vol = %f\n",vol);
+			rad = 0.02 * vol/pow(cg[0].npat, 1.0/3.0);
+//printf("~1 rad = %f\n",rad);
+
+			if (dovrml == 3)	// Hack
+				rad = 0.02;
+			else if (dovrml == 4)	// Hack
+				rad = 0.015;
+		}
+
+		if (dovrml && (dovrml == 3 || dovrml == 4)) {		/* RGB/YCC device plot */
+			if (!cg[0].isrgb || !cg[1].isrgb)
+				error("Both files must have RGB devices space for -d option");
 		}
 
 		/* Do overall results */
@@ -932,13 +1185,45 @@ int main(int argc, char *argv[])
 				merr = de;
 
 			if (dovrml) {
-				if (de > 1e-6) {
-					wrl->add_vertex(wrl, 0, cg[0].pat[j].v);
-					wrl->add_vertex(wrl, 0, cg[1].pat[j].v);
-				}
-				if (dovrml == 2) {
-					wrl->add_marker(wrl, cg[0].pat[j].v, NULL, rad);
-					wrl->add_marker(wrl, cg[1].pat[j].v, NULL, rad);
+				if ((dovrml == 3 || dovrml == 4)) {		/* RGB/YCC device plot */
+					double *val1, *val2;
+					int k;
+
+					if (dovrml == 3) {
+						val1 = cg[0].pat[i].rgb; 
+						val2 = cg[1].pat[match[i]].rgb;
+					} else {
+						val1 = cg[0].pat[i].ycc; 
+						val2 = cg[1].pat[match[i]].ycc;
+					}
+
+					de = icmNorm33(val1, val2);
+
+					if (de > 1e-6) {
+						wrl->add_vertex(wrl, 0, val1);
+						wrl->add_vertex(wrl, 0, val2);
+					}
+
+#ifdef NEVER	// Green target
+					wrl->add_marker(wrl, val1, green, rad);
+
+#else		// Natural color
+					for (k = 0; k < 3; k++)
+						col[k] = 0.3 + 0.7 * (cg[0].pat[i].rgb[k] - min[k])/(max[k] - min[k]);
+					wrl->add_marker(wrl, val1, col, rad);
+#endif
+
+					wrl->add_marker_trans(wrl, val2, red, 0.3, rad * 0.99);
+
+				} else {		/* PCS */
+					if (de > 1e-6) {
+						wrl->add_vertex(wrl, 0, cg[0].pat[i].v);
+						wrl->add_vertex(wrl, 0, cg[1].pat[match[i]].v);
+					}
+					if (dovrml == 2) {
+						wrl->add_marker(wrl, cg[0].pat[i].v, green, rad);
+						wrl->add_marker(wrl, cg[1].pat[match[i]].v, red, rad);
+					}
 				}
 			}
 
@@ -997,6 +1282,18 @@ int main(int argc, char *argv[])
 			fprintf(verbo,"No of test patches in best 90%% are = %d\n",n90);
 		}
 		printf("Verify results:\n");
+		if (norm == 4)
+			printf("  L*a*b* ref. = average XYZ %f %f %f\n",cg[0].w[0],cg[0].w[1],cg[0].w[2]);
+		else if (norm == 1) {
+			printf("  File 1 L* ref. Y %f\n", cg[0].w[1]);
+			printf("  File 2 L* ref. Y %f\n", cg[1].w[1]);
+		} else if (norm == 2) {
+			printf("  File 1 L*a*b* ref. XYZ %f %f %f\n", cg[0].w[0],cg[0].w[1],cg[0].w[2]);
+			printf("  File 2 L*a*b* ref. XYZ %f %f %f\n", cg[1].w[0],cg[1].w[1],cg[1].w[2]);
+		} else if (norm == 3) {
+			printf("  File 1 L* ref. X+Y+Z %f %f %f\n", cg[0].w[0],cg[0].w[1],cg[0].w[2]);
+			printf("  File 2 L* ref. X+Y+Z %f %f %f\n", cg[1].w[0],cg[1].w[1],cg[1].w[2]);
+		}
 		printf("  Total errors%s:     peak = %f, avg = %f\n", cie2k ? " (CIEDE2000)" : cie94 ? " (CIE94)" : "", merr, aerr);
 		printf("  Worst 10%% errors%s: peak = %f, avg = %f\n", cie2k ? " (CIEDE2000)" : cie94 ? " (CIE94)" : "", merr10, aerr10);
 		printf("  Best  90%% errors%s: peak = %f, avg = %f\n", cie2k ? " (CIEDE2000)" : cie94 ? " (CIE94)" : "", merr90, aerr90);
diff --git a/profile/invprofcheck.c b/profile/invprofcheck.c
index 0965096..8cb9f34 100644
--- a/profile/invprofcheck.c
+++ b/profile/invprofcheck.c
@@ -37,6 +37,8 @@
 #include "numlib.h"
 #include "icc.h"
 #include "xicc.h"
+#include "vrml.h"
+#include "ui.h"
 
 /* Resolution of the sampling modes */
 #define TRES 11
@@ -92,19 +94,14 @@ void usage(void) {
 	fprintf(stderr," -R res       Specific grid resolution\n");
 	fprintf(stderr," -c           Show CIE94 delta E values\n");
 	fprintf(stderr," -k           Show CIEDE2000 delta E values\n");
-	fprintf(stderr," -w           create VRML visualisation (profile.wrl)\n");
-	fprintf(stderr," -x           Use VRML axes\n");
+	fprintf(stderr," -w           create %s visualisation (profile%s)\n",vrml_format(),vrml_ext());
+	fprintf(stderr," -x           Use %s axes\n",vrml_format());
 	fprintf(stderr," -e           Color vectors acording to delta E\n");
 	fprintf(stderr," profile.icm  Profile to check\n");
 	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 make_de_lines(FILE *wrl);
-void end_vrml(FILE *wrl);
+static void DE2RGB(double *out, double in);
 
 #if defined(__IBMC__) && defined(_M_IX86)
 void bug_workaround(int *co) { };			/* Workaround optimiser bug */
@@ -123,14 +120,14 @@ main(
 	int doaxes = 0;
 	int dodecol = 0;
 	char in_name[MAXNAMEL+1];
-	char out_name[MAXNAMEL+1], *xl;		/* VRML name */
+	char out_name[MAXNAMEL+1], *xl;		/* VRML/X3D name */
 	icmFile *rd_fp;
 	icc *icco;
 	int rv = 0;
 	int tres = TRES;
 	double tlimit = -1.0;
 	double klimit = -1.0;
-	FILE *wrl = NULL;
+	vrml *wrl = NULL;
 
 	error_program = "invprofcheck";
 
@@ -206,7 +203,7 @@ main(
 				klimit = limit/100.0;
 			}
 
-			/* VRML */
+			/* VRML/X3D */
 			else if (argv[fa][1] == 'w' || argv[fa][1] == 'W')
 				dovrml = 1;
 
@@ -242,7 +239,7 @@ main(
 	strncpy(out_name,in_name,MAXNAMEL-4); out_name[MAXNAMEL-4] = '\000';
 	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)
@@ -282,8 +279,8 @@ main(
 		}
 
 		if (dovrml) {
-			wrl = start_vrml(out_name, doaxes);
-			start_line_set(wrl);
+			wrl = new_vrml(out_name, doaxes, vrml_lab);
+			wrl->start_line_set(wrl, 0);
 		}
 		
 		/* Grab any device calibration curves */
@@ -357,8 +354,21 @@ main(
 
 				/* Delta E */
 				if (dovrml) {
-					add_vertex(wrl, pcsin);
-					add_vertex(wrl, pcsout);
+					int ix[2];
+
+					/* Add the verticies */
+					ix[0] = wrl->add_vertex(wrl, 0, pcsin);
+					ix[1] = wrl->add_vertex(wrl, 0, pcsout);
+
+					/* Add the line */
+					if (dodecol) {		/* Lines with color determined by length */
+						double rgb[3];
+						DE2RGB(rgb, icmNorm33(pcsin, pcsout));
+						wrl->add_col_line(wrl, 0, ix, rgb);
+	
+					} else {	/* Natural color */
+						wrl->add_line(wrl, 0, ix);
+					}
 				}
 	
 				/* Check the result */
@@ -386,11 +396,8 @@ main(
 			}
 		}
 		if (dovrml) {
-			if (dodecol)
-				make_de_lines(wrl);
-			else
-				make_lines(wrl, 2);
-			end_vrml(wrl);
+			wrl->make_lines_vc(wrl, 0, 0.0);
+			wrl->del(wrl);
 		}
 
 		printf("Profile check complete, errors%s: max. = %f, avg. = %f, RMS = %f\n",
@@ -409,291 +416,8 @@ main(
 
 
 /* ------------------------------------------------ */
-/* Some simple functions to do basix VRML work */
-/* !!! Should change to plot/vrml lib !!! */
-
-#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);
-static void DE2RGB(double *out, double in);
-
-FILE *start_vrml(char *name, int doaxes) {
-	FILE *wrl;
-
-	/* Define the axis boxes */
-	struct {
-		double x, y, z;			/* Box center */
-		double wx, wy, wz;		/* Box size */
-		double r, g, b;			/* Box color */
-	} 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 */
-	};
-
-	/* Define the labels */
-	struct {
-		double x, y, z;
-		double size;
-		char *string;
-		double r, g, b;
-	} labels[6] = {
-		{ -2, 2, -GAMUT_LCENT + 100 + 10, 10, "+L*",  .7, .7, .7 },	/* Top of L axis */
-		{ -2, 2, -GAMUT_LCENT - 10,      10, "0",    .7, .7, .7 },	/* Bottom of L axis */
-		{ 100 + 5, -3,  0-GAMUT_LCENT,  10, "+a*",  1,  0,  0 },	/* +a (red) axis */
-		{ -5, -100 - 10, 0-GAMUT_LCENT,  10, "-b*",  0,  0,  1 },	/* -b (blue) axis */
-		{ -100 - 15, -3, 0-GAMUT_LCENT,  10, "-a*",  0,  0,  1 },	/* -a (green) axis */
-		{ -5,  100 + 5, 0-GAMUT_LCENT,  10, "+b*",  1,  1,  0 },	/* +b (yellow) axis */
-	};
-
-	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) {
-		int n;
-		fprintf(wrl,"    # Lab axes as boxes:\n");
-		for (n = 0; n < 5; n++) {
-			fprintf(wrl,"    Transform { translation %f %f %f\n", axes[n].x, axes[n].y, axes[n].z);
-			fprintf(wrl,"      children [\n");
-			fprintf(wrl,"        Shape{\n");
-			fprintf(wrl,"          geometry Box { size %f %f %f }\n",
-			                       axes[n].wx, axes[n].wy, axes[n].wz);
-			fprintf(wrl,"          appearance Appearance { material Material ");
-			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[n].r, axes[n].g, axes[n].b);
-			fprintf(wrl,"        }\n");
-			fprintf(wrl,"      ]\n");
-			fprintf(wrl,"    }\n");
-		}
-		fprintf(wrl,"    # Axes identification:\n");
-		for (n = 0; n < 6; n++) {
-			fprintf(wrl,"    Transform { translation %f %f %f\n", labels[n].x, labels[n].y, labels[n].z);
-			fprintf(wrl,"      children [\n");
-			fprintf(wrl,"        Shape{\n");
-			fprintf(wrl,"          geometry Text { string [\"%s\"]\n",labels[n].string);
-			fprintf(wrl,"            fontStyle FontStyle { family \"SANS\" style \"BOLD\" size %f }\n",
-			                                  labels[n].size);
-			fprintf(wrl,"                        }\n");
-			fprintf(wrl,"          appearance Appearance { material Material ");
-			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", labels[n].r, labels[n].g, labels[n].b);
-			fprintf(wrl,"        }\n");
-			fprintf(wrl,"      ]\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");
-}
-
-/* Assume 2 ppset, and make line color prop to length */
-void make_de_lines(FILE *wrl) {
-	int i, j;
-
-	fprintf(wrl,"      ]\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"  coordIndex [\n");
-
-	for (i = 0; i < npoints;) {
-		for (j = 0; j < 2; 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], ss;
-		for (ss = 0.0, j = 0; j < 3; j++) {
-			double tt = (pary[i & ~1].pp[j] - pary[i | 1].pp[j]);
-			ss += tt * tt;
-		}
-		ss = sqrt(ss);
-		DE2RGB(rgb, ss);
-		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;
-}
-
 /* Convert a delta E value into a signal color: */
-static void
-DE2RGB(double *out, double in) {
+static void DE2RGB(double *out, double in) {
 	struct {
 		double de;
 		double r, g, b;
diff --git a/profile/ls2ti3.c b/profile/ls2ti3.c
new file mode 100644
index 0000000..3330e16
--- /dev/null
+++ b/profile/ls2ti3.c
@@ -0,0 +1,380 @@
+
+/* 
+ * Argyll Color Correction System
+ *
+ * Read in the RGB & CIE data from a LightSpace format .bcs XML file
+ * and convert it into a .ti3 CGATs format suitable for the Argyll CMS.
+ *
+ * Derived from  txt2cgats.c 
+ * Author: Graeme W. Gill
+ * Date:   16/11/00
+ *
+ * Copyright 2000 - 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.
+ */
+
+/* TTBD
+	
+ */
+
+#define DEBUG
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <sys/types.h>
+#include <time.h>
+#include <string.h>
+#include <stdarg.h>
+#include "copyright.h"
+#include "aconfig.h"
+#include "cgats.h"
+#include "xspect.h"
+#include "insttypes.h"
+#include "mxml.h"
+#include "numlib.h"
+#include "ui.h"
+
+#define DEB 6
+
+void
+usage(char *mes) {
+	fprintf(stderr,"Convert LightSpace raw RGB device profile data to Argyll CGATS data, Version %s\n",ARGYLL_VERSION_STR);
+	fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
+	if (mes != NULL)
+		fprintf(stderr,"error: %s\n",mes);
+	fprintf(stderr,"usage: ls2ti3 [-v] infile outbase\n");
+/*	fprintf(stderr," -v            Verbose mode\n"); */
+	fprintf(stderr," infile        Input LightSpace .bcs file\n");
+	fprintf(stderr," outbasename   Output file basename for .ti3\n");
+	exit(1);
+	}
+
+/* XML data type callback for mxmlLoadFile() */
+mxml_type_t
+type_cb(mxml_node_t *node) {
+	mxml_node_t *parent = mxmlGetParent(node);
+	const char *pname;
+	const char *name = node->value.element.name;
+
+	if (parent == NULL)
+		return MXML_TEXT;
+
+	pname = parent->value.element.name;
+
+//	printf("~1 type_cb got pnode '%s' node '%s'\n",pname, name);
+
+	if (strcmp(pname, "stimuli") == 0
+	 && (strcmp(name, "red") == 0
+	  || strcmp(name, "green") == 0
+	  || strcmp(name, "blue") == 0))
+		return MXML_REAL;
+
+	if (strcmp(pname, "XYZ") == 0
+	 && (strcmp(name, "X") == 0
+	  || strcmp(name, "Y") == 0
+	  || strcmp(name, "Z") == 0))
+		return MXML_REAL;
+
+#ifdef NEVER
+	if (strcmp(pname, "FileInformation") == 0
+	 && (strcmp(name, "Creator") == 0
+	  || strcmp(name, "CreationDate") == 0
+	  || strcmp(name, "Description") == 0))
+		return MXML_OPAQUE;		/* Don't split strings up */
+#endif
+
+	return MXML_TEXT;
+}
+
+struct _patch {
+	int pno;				/* Patch number */
+
+	double dev[MAX_CHAN];	/* Device value */
+	double XYZ[3];			/* XYZ value */
+
+}; typedef struct _patch patch;
+
+
+int main(int argc, char *argv[]) {
+	int i, j;
+	int fa,nfa;				/* current argument we're looking at */
+	int verb = 0;
+	static char inname[MAXNAMEL+1] = { 0 };		/* Input LightSpace .xml file */
+	static char outname[MAXNAMEL+1] = { 0 };		/* Output cgats .ti3 file base name */
+	double dev_scale = 1.0;	/* Device value scaling */
+
+    FILE *ifp;
+    mxml_node_t *tree, *pnode, *node;
+    mxml_node_t *top, *data;
+	const char *attr, *name;
+	patch *patches;
+
+	cgats *ocg;				/* output cgats structure for .ti3 */
+	time_t clk = time(0);
+	struct tm *tsp = localtime(&clk);
+	char *atm = asctime(tsp); /* Ascii time */
+
+	int islab = 0;			/* CIE is Lab rather than XYZ */
+
+	int npat = 0;			/* Number of patches */
+	int ndchan = 3;			/* Number of device channels, 0 = no device, RGB = 3, CMYK = 4 */
+	cgats_set_elem *setel;	/* Array of set value elements */
+
+	error_program = "ls2ti3";
+
+	if (argc < 3)
+		usage("Too few arguments");
+
+	/* Process the arguments */
+	for(fa = 1;fa < argc;fa++) {
+		nfa = fa;					/* skip to nfa if next argument is used */
+		if (argv[fa][0] == '-') {	/* Look for any flags */
+			char *na = NULL;		/* next argument after flag, null if none */
+
+			if (argv[fa][2] != '\000')
+				na = &argv[fa][2];		/* next is directly after flag */
+			else {
+				if ((fa+1) < argc) {
+					if (argv[fa+1][0] != '-') {
+						nfa = fa + 1;
+						na = argv[nfa];		/* next is seperate non-flag argument */
+					}
+				}
+			}
+
+			if (argv[fa][1] == '?')
+				usage(NULL);
+
+			else if (argv[fa][1] == 'v' || argv[fa][1] == 'V')
+				verb = 1;
+			else 
+				usage("Unknown flag");
+		} else
+			break;
+	}
+
+	/* Get the file name argument */
+	if (fa >= argc || argv[fa][0] == '-') usage("input filename not found");
+	strncpy(inname,argv[fa++],MAXNAMEL); inname[MAXNAMEL] = '\000';
+
+	if (fa >= argc || argv[fa][0] == '-') usage("output basename not found");
+	strncpy(outname,argv[fa++],MAXNAMEL-4); outname[MAXNAMEL-4] = '\000';
+	strcat(outname,".ti3");
+
+    if ((ifp = fopen(inname, "r")) == NULL)
+		error("Unable to read '%s'",inname);
+
+#ifndef NEVER
+    tree = mxmlLoadFile(NULL, ifp, type_cb);
+#else
+	tree = mxmlSAXLoadFd(NULL, fileno(ifp), type_cb, sax_cb, (void *)p);
+#endif
+    fclose(ifp);
+
+	if (tree == NULL)
+		error("Parsing '%s' failed",inname);
+
+	if ((top = mxmlFindElement(tree, tree, NULL, NULL, NULL, MXML_DESCEND_FIRST)) == NULL
+	 || mxmlGetType(top) != MXML_ELEMENT)
+		error("Failed to find top element in '%s'",inname);
+
+	if (strcmp(top->value.element.name, "builder_color_space") != 0)
+		error("'%s' doesn't seem to be a LightSpace .bcs file ?",inname);
+
+#ifdef NEVER
+	/* Check that its CxF3 */
+	if ((attr = mxmlElementGetAttr(top, "xmlns")) == NULL)
+		error("Failed to find CxF attribute %s in '%s'","xmlns", inname);
+
+	if (strcmp(attr, "http://colorexchangeformat.com/CxF3-core") != 0)
+		error("File '%s' is not CxF format",inname);
+
+	/* Look for header information */
+	if ((pnode = mxmlFindElement(top, top, "head", NULL, NULL, MXML_DESCEND_FIRST)) == NULL)
+		error("Failed to find head in '%s'",inname);
+
+	/* Grab the creation date */
+	if ((node = mxmlFindElement(pnode, pnode, "created", NULL, NULL, MXML_DESCEND_FIRST)) == NULL)
+		error();
+
+	name = mxmlGetOpaque(node);
+#endif
+
+	/* Locate the data node */
+	if ((data = mxmlFindElement(top, top, "data", NULL, NULL, MXML_DESCEND_FIRST)) == NULL)
+		error("Failed to find data in '%s'",inname);
+
+	/* Get the number of patches */
+	if ((attr = mxmlElementGetAttr(data, "frames")) == NULL)
+		error("Failed to find data attribute 'frames' in '%s'", inname);
+
+	npat = atoi(attr);
+
+	if (npat <= 0)
+		error("Illegal number of patches %d in '%s",inname);
+
+	a1logd(g_log, DEB, "npat = %d\n",npat);
+
+	if ((patches = calloc(npat, sizeof(patch))) == NULL)
+		error("malloc failed");
+
+	/* Read all the patches */
+	node = mxmlFindElement(data, data, "patch", NULL, NULL, MXML_DESCEND_FIRST);
+	for (i = 0; node != NULL && i < npat;) {
+		int j;
+	    mxml_node_t *stim, *res, *xyz, *val;
+		char *rgb_key[3] = { "red", "green", "blue" };
+		char *xyz_key[3] = { "X", "Y", "Z" };
+		
+		if (mxmlGetType(node) != MXML_ELEMENT) {
+			a1logd(g_log, DEB, "skipping non element node type %d\n",mxmlGetType(node));
+			goto next;
+		}
+		a1logd(g_log, DEB, "read node '%s'\n",node->value.element.name);
+
+		if ((attr = mxmlElementGetAttr(node, "frame")) == NULL) {
+			a1logd(g_log, DEB, "read_cxf: skipping node without frame\n");
+			goto next;							/* Skip this one */
+		}
+		patches[i].pno = atoi(attr);
+
+		a1logd(g_log, DEB, "got patch %d no %d\n",i,patches[i].pno);
+
+		/* Read the RGB stimulus */
+		if ((stim = mxmlFindElement(node, node, "stimuli", NULL, NULL, MXML_DESCEND_FIRST))
+			                                                                       == NULL) {
+			a1logd(g_log, DEB, "Can't find 'stimuli' in patch %d - skipping\n",patches[i].pno);
+			goto next;
+		}
+
+		for (j = 0; j < 3; j++) {
+			if ((val = mxmlFindElement(stim, stim, rgb_key[j], NULL, NULL, MXML_DESCEND_FIRST))
+				                                                                       == NULL) {
+				a1logd(g_log, DEB, "failed to find RGB component %s\n",rgb_key[j]);
+				goto next;
+			}
+			patches[i].dev[j] = mxmlGetReal(val);
+			a1logd(g_log, DEB, "got RGB component %s value %f\n",rgb_key[j],patches[i].dev[j]);
+		}
+
+		/* Read the XYZ results */
+		if ((res = mxmlFindElement(node, node, "results", NULL, NULL, MXML_DESCEND_FIRST))
+			                                                                       == NULL) {
+			a1logd(g_log, DEB, "Can't find 'results' in patch %d - skipping\n",patches[i].pno);
+			goto next;
+		}
+
+		if ((xyz = mxmlFindElement(res, res, "XYZ", NULL, NULL, MXML_DESCEND_FIRST))
+			                                                                       == NULL) {
+			a1logd(g_log, DEB, "Can't find 'XYZ' in patch %d - skipping\n",patches[i].pno);
+			goto next;
+		}
+
+		for (j = 0; j < 3; j++) {
+			if ((val = mxmlFindElement(xyz, xyz, xyz_key[j], NULL, NULL, MXML_DESCEND_FIRST))
+				                                                                       == NULL) {
+				a1logd(g_log, DEB, "failed to find XYZ component %s\n",xyz_key[j]);
+				goto next;
+			}
+			patches[i].XYZ[j] = mxmlGetReal(val);
+			a1logd(g_log, DEB, "got XYZ component %s value %f\n",xyz_key[j],patches[i].XYZ[j]);
+		}
+		i++;
+
+		/* Add patch to output CGATS */
+
+	next:;	/* Next color */
+		node = mxmlGetNextSibling(node);
+	}
+
+	mxmlDelete(tree);
+
+
+	/* Setup output cgats file */
+	ocg = new_cgats();	/* Create a CGATS structure */
+	ocg->add_other(ocg, "CTI3"); 	/* our special type is Calibration Target Information 3 */
+	ocg->add_table(ocg, tt_other, 0);	/* Start the first table */
+
+	ocg->add_kword(ocg, 0, "DESCRIPTOR", "Argyll Calibration Target chart information 3",NULL);
+	ocg->add_kword(ocg, 0, "ORIGINATOR", "Argyll target", NULL);
+	atm[strlen(atm)-1] = '\000';	/* Remove \n from end */
+	ocg->add_kword(ocg, 0, "CREATED",atm, NULL);
+	ocg->add_kword(ocg, 0, "DEVICE_CLASS","DISPLAY", NULL);	/* What sort of device this is */
+
+	/* Note what instrument the chart was read with */
+	/* Assume this - could try reading from file INSTRUMENTATION "SpectroScan" ?? */
+	ocg->add_kword(ocg, 0, "TARGET_INSTRUMENT", inst_name(instSpectrolino) , NULL);
+
+	/* Don't make any assumptions about normalisation (this is default anyway) */
+	ocg->add_kword(ocg, 0, "NORMALIZED_TO_Y_100","NO", NULL);
+
+	/* Fields we want */
+	ocg->add_field(ocg, 0, "SAMPLE_ID", nqcs_t);
+
+	if (ndchan == 3) {
+		ocg->add_field(ocg, 0, "RGB_R", r_t);
+		ocg->add_field(ocg, 0, "RGB_G", r_t);
+		ocg->add_field(ocg, 0, "RGB_B", r_t);
+		if (islab)
+			ocg->add_kword(ocg, 0, "COLOR_REP","RGB_LAB", NULL);
+		else
+			ocg->add_kword(ocg, 0, "COLOR_REP","RGB_XYZ", NULL);
+	}
+
+	if (islab) {
+		ocg->add_field(ocg, 0, "LAB_L", r_t);
+		ocg->add_field(ocg, 0, "LAB_A", r_t);
+		ocg->add_field(ocg, 0, "LAB_B", r_t);
+	} else {
+		ocg->add_field(ocg, 0, "XYZ_X", r_t);
+		ocg->add_field(ocg, 0, "XYZ_Y", r_t);
+		ocg->add_field(ocg, 0, "XYZ_Z", r_t);
+	}
+
+	/* LS uses 0.0 .. 1.0 */
+	dev_scale = 100.0/1.0;
+
+	/* Write out the patch info to the output CGATS file */
+	if ((setel = (cgats_set_elem *)malloc(
+	     sizeof(cgats_set_elem) * ocg->t[0].nfields)) == NULL)
+		error("Malloc failed!");
+
+	/* Write out the patch info to the output CGATS file */
+	for (i = 0; i < npat; i++) {
+		char id[100];
+		int k = 0;
+
+		/* SAMPLE ID */
+		sprintf(id, "%d", patches[i].pno);
+		setel[k++].c = id;
+
+		if (ndchan == 3) {
+			setel[k++].d = dev_scale * patches[i].dev[0];
+			setel[k++].d = dev_scale * patches[i].dev[1];
+			setel[k++].d = dev_scale * patches[i].dev[2];
+		} 
+
+		setel[k++].d = patches[i].XYZ[0];
+		setel[k++].d = patches[i].XYZ[1];
+		setel[k++].d = patches[i].XYZ[2];
+
+		ocg->add_setarr(ocg, 0, setel);
+	}
+
+	free(setel);
+
+	if (ocg->write_name(ocg, outname))
+		error("Write error : %s",ocg->err);
+
+	/* Clean up */
+	ocg->del(ocg);
+
+	return 0;
+}
+
+
+
diff --git a/profile/mppcheck.c b/profile/mppcheck.c
index 2efe878..0fcfbb7 100644
--- a/profile/mppcheck.c
+++ b/profile/mppcheck.c
@@ -33,6 +33,7 @@
 #include "xicc.h"
 #include "prof.h"
 #include "sort.h"
+#include "ui.h"
 
 void
 usage(void) {
diff --git a/profile/mppprof.c b/profile/mppprof.c
index 0c3593f..74d3a71 100644
--- a/profile/mppprof.c
+++ b/profile/mppprof.c
@@ -62,7 +62,8 @@
 #include "numlib.h"
 #include "xicc.h"
 #include "prof.h"
-#include "../h/sort.h"
+#include "sort.h"
+#include "ui.h"
 
 void
 usage(void) {
diff --git a/profile/printcal.c b/profile/printcal.c
index f9107eb..38ceb9b 100644
--- a/profile/printcal.c
+++ b/profile/printcal.c
@@ -51,6 +51,7 @@
 #include "rspl.h"
 #include "xicc.h"
 #include "plot.h"
+#include "ui.h"
 
 
 #define RSPLFLAGS (0 /* | RSPL_2PASSSMTH | RSPL_EXTRAFIT2 */)
@@ -1366,7 +1367,22 @@ int main(int argc, char *argv[]) {
 		           NULL);				/* iwidth */
 
 
-		if (verb > 1) {
+		/* Compute & show fit quality */
+		if (verb > 0) {
+			double avgde = 0.0, maxde = 0.0;
+			for (i = 0; i < n_pvals[j]; i++) {
+				co tp;	/* Test point */
+				double de;
+				tp.p[0] = pvals[j][i].dev;
+				raw[j]->interp(raw[j], &tp);
+				de = icmLabDE(pvals[j][i].Lab, tp.v);
+
+				avgde += de;
+				if (de > maxde)
+					maxde = de;
+			}
+			avgde /= (double)n_pvals[j];
+			printf("Chan %d raw fit avg DE %f, max %f\n",j,avgde,maxde);
 		}
 
 		free(dpoints);
diff --git a/profile/prof.h b/profile/prof.h
index 6013d99..098c03f 100644
--- a/profile/prof.h
+++ b/profile/prof.h
@@ -50,6 +50,7 @@ void make_output_icc(
 	int verify,				/* nz to print verification */
 	int clipprims,			/* Clip white, black and primaries */
 	double wpscale,			/* >= 0.0 for media white point scale factor */
+//	double *bpo,			/* != NULL for XYZ black point override */
 	icxInk *ink,			/* Ink limit/black generation setup */
 	char *in_name,			/* input .ti3 file name */
 	char *file_name,		/* output icc name */
diff --git a/profile/profcheck.c b/profile/profcheck.c
index 1895713..32cddad 100644
--- a/profile/profcheck.c
+++ b/profile/profcheck.c
@@ -27,6 +27,8 @@
 
 #undef DEBUG
 
+#undef HACK					/* Print per patch XYZ differences */
+
 #define IMP_MONO			/* Turn on development code */
 
 #define verbo stdout
@@ -42,6 +44,9 @@
 #include "xicc.h"
 #include "insttypes.h"
 #include "sort.h"
+#include "plot.h"
+#include "vrml.h"
+#include "ui.h"
 
 void
 usage(void) {
@@ -51,10 +56,13 @@ usage(void) {
 	fprintf(stderr," -v [level]      Verbosity level (default 1), 2 to print each DE\n");
 	fprintf(stderr," -c              Show CIE94 delta E values\n");
 	fprintf(stderr," -k              Show CIEDE2000 delta E values\n");
-	fprintf(stderr," -w              create VRML visualisation (iccprofile.wrl)\n");
-	fprintf(stderr," -x              Use VRML axes\n");
-	fprintf(stderr," -m              Make VRML lines a minimum of 0.5\n");
+	fprintf(stderr," -w              create %s visualisation (iccprofile%s)\n",vrml_format(),vrml_ext());
+	fprintf(stderr," -x              Use %s axes\n",vrml_format());
+	fprintf(stderr," -m              Make %s lines a minimum of 0.5\n",vrml_format());
 	fprintf(stderr," -e              Color vectors acording to delta E\n");
+	fprintf(stderr," -h              Plot a histogram of delta E's\n");
+	fprintf(stderr," -s              Sort output by delta E\n");
+	fprintf(stderr," -P N.NN         Create a pruned .ti3 with points less or equal to N.NN delta E\n");
 	fprintf(stderr," -d devval1,deval2,devvalN\n");
 	fprintf(stderr,"                 Specify a device value to sort against\n");
 	fprintf(stderr," -p              Sort device value by PCS (Lab) target\n");
@@ -70,12 +78,7 @@ usage(void) {
 	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 make_de_lines(FILE *wrl);
-void end_vrml(FILE *wrl);
+static void DE2RGB(double *out, double in);
 
 /* Patch value type */
 typedef struct {
@@ -83,10 +86,21 @@ typedef struct {
 	char slo[50];		/* sample location, "" if not known */
 	double p[MAX_CHAN];	/* Device value */
 	double v[3];		/* CIE value */
-	double dp;			/* Delta from target value */
-	double dv;			/* Delta from CIE value */
+
+	double pv[3];		/* Profile CIE value */
+	double de;			/* Delta E to profile CIE value */
+
+	double dp;			/* Delta p[] from target device value */
+	double dv;			/* Delta E from CIE value */
 } pval;
 
+/* Histogram bin type */
+typedef struct {
+	int count;			/* Raw count */
+	double val;			/* Normalized value */
+	double min, max;	/* Bin range */
+} hbin;
+
 int main(int argc, char *argv[])
 {
 	int fa,nfa;				/* current argument we're looking at */
@@ -104,8 +118,8 @@ int main(int argc, char *argv[])
 	icRenderingIntent intent = icAbsoluteColorimetric;
 	icc *rd_icco;
 	icmLuBase *luo;
-	char out_name[MAXNAMEL+1], *xl;		/* VRML name */
-	FILE *wrl = NULL;
+	char out_name[MAXNAMEL+1], *xl;		/* VRML/X3D name */
+	vrml *wrl = NULL;
 
 	int fwacomp = 0;			/* FWA compensation */
 	int isdisp = 0;				/* nz if this is a display device, 0 if output */
@@ -117,12 +131,18 @@ int main(int argc, char *argv[])
 	xspect cust_illum;			/* Custom illumination spectrum */
 	icxObserverType observ = icxOT_CIE_1931_2;
 
+	int sortbyde = 0;			/* Sort by delta E */
+
 	int ddevv = 0;				/* Do device value sort */
 	double devval[MAX_CHAN];	/* device value to sort on */
-	int sortbypcs = 0;			/* Sort by PCS */
+	int sortbypcs = 0;			/* Sort by PCS error to device target */
+
+	int dohisto = 0;			/* Plot histogram of delta E's */
+	double prune = 0.0;			/* If > 0.0, created a pruned .ti3 file */
 
 	int npat;					/* Number of patches */
 	pval *tpat;					/* Patch input values */
+	pval **stpat;				/* Pointers to internal sorted tpat[] */
 	int i, j, rv = 0;
 	icColorSpaceSignature devspace = 0;	/* The device colorspace */
 	int isAdditive = 0;			/* 0 if subtractive, 1 if additive colorspace */
@@ -168,7 +188,7 @@ int main(int argc, char *argv[])
 				}
 			}
 
-			/* VRML */
+			/* VRML/X3D */
 			else if (argv[fa][1] == 'w')
 				dovrml = 1;
 
@@ -194,6 +214,23 @@ int main(int argc, char *argv[])
 				cie2k = 1;
 			}
 
+			/* Plot histogram */
+			else if (argv[fa][1] == 'h')
+				dohisto = 1;
+
+			/* Sort by delta E */
+			else if (argv[fa][1] == 's')
+				sortbyde = 1;
+
+			/* Create a pruned .ti3 */
+			else if (argv[fa][1] == 'P') {
+				fa = nfa;
+				if (na == NULL) usage();
+				prune = atof(na);
+				if (prune <= 0.0)
+					usage();
+			}
+
 			/* Device sort value */
 			else if (argv[fa][1] == 'd') {
 				char *tp, buf[200];
@@ -354,7 +391,7 @@ int main(int argc, char *argv[])
 	strncpy(out_name,iccname,MAXNAMEL-4); out_name[MAXNAMEL-4] = '\000';
 	if ((xl = strrchr(out_name, '.')) == NULL)	/* Figure where extention is */
 		xl = out_name + strlen(out_name);
-	strcpy(xl,".wrl");
+	xl[0] = '\000';			/* Remove extension */
 
 	if (fwacomp && spec == 0)
 		error ("FWA compensation only works when viewer and/or illuminant selected");
@@ -504,6 +541,12 @@ int main(int argc, char *argv[])
 	if ((tpat = (pval *)malloc(sizeof(pval) * npat)) == NULL)
 		error("Malloc failed - tpat[]");
 
+	if ((stpat = (pval **)malloc(sizeof(pval *) * npat)) == NULL)
+		error("Malloc failed - stpat[]");
+
+	for (i = 0; i < npat; i++)
+		stpat[i] = &tpat[i];
+
 	/* Read in the CGATs fields */
 	{
 		int sidx;					/* Sample ID index */
@@ -840,11 +883,212 @@ int main(int argc, char *argv[])
 			}
 		}
 
-		icg->del(icg);		/* Clean up */
 	}	/* End of reading in CGATs file */
 
 	/* - - - - - - - - - - */
-	/* Check the forward profile accuracy against the data points */
+	/* Compute the delta E of each point against the profile value */
+
+	/* Open up the file for reading */
+	if ((rd_fp = new_icmFileStd_name(iccname,"r")) == NULL)
+		error("Write: Can't open file '%s'",iccname);
+
+	if ((rd_icco = new_icc()) == NULL)
+		error("Read: Creation of ICC object failed");
+
+	/* Read the header and tag list */
+	if ((rv = rd_icco->read(rd_icco,rd_fp,0)) != 0)
+		error("Read: %d, %s",rv,rd_icco->err);
+
+	/* Get the Fwd table, absolute with Lab override */
+	if ((luo = rd_icco->get_luobj(rd_icco, icmFwd, intent,
+	                              icSigLabData, icmLuOrdNorm)) == NULL) {
+		error("%d, %s",rd_icco->errc, rd_icco->err);
+	}
+
+	for (i = 0; i < npat; i++) {
+
+		/* Lookup the patch value in the profile */
+		if (luo->lookup(luo, tpat[i].pv, tpat[i].p) > 1)
+			error("%d, %s",rd_icco->errc,rd_icco->err);
+
+		if (cie2k)
+			tpat[i].de = icmCIE2K(tpat[i].v, tpat[i].pv);
+		else if (cie94)
+			tpat[i].de = icmCIE94(tpat[i].v, tpat[i].pv);
+		else
+			tpat[i].de = icmLabDE(tpat[i].v, tpat[i].pv);
+	}
+
+	/* - - - - - - - - - - */
+	/* Sort by delta E */
+	if (sortbyde) {
+		/* Sort the dE's */
+#define HEAP_COMPARE(A,B) (A.de > B.de)
+		HEAPSORT(pval, tpat, npat);
+#undef HEAP_COMPARE
+	}
+
+	/* - - - - - - - - - - */
+	/* Plot a dE histogram */
+	if (dohisto) {
+		double demax = -1e6, demin = 1e6;
+		int maxbins = 50;		/* Maximum bins */
+		int minbins = 20;		/* Target minimum bins (depends on distribution) */ 
+		int mincount = 10;		/* Minimum number of points in a bin */
+		double mbwidth;
+		int nbins = 0;
+		hbin *bins = NULL;
+		double tval;
+		double *x, *y;
+		
+		/* Figure out the range of dE's */
+		for (i = 0; i < npat; i++) {
+			double de = tpat[i].de;
+
+			if (de > demax)
+				demax = de;
+			if (de < demin)
+				demin = de;
+		}
+
+		if (demax < 1e-6)
+			error("histogram: dE range is too small to plot");
+
+		/* Bin width that gives maxbins */
+		mbwidth = demax / maxbins;
+		
+		/* Reduce mincount if needed to get minbins */
+		if (npat/minbins < mincount)
+			mincount = npat/minbins;
+
+		if ((bins = (hbin *)calloc(maxbins, sizeof(hbin))) == NULL)
+			error("malloc of histogram bins failed");
+
+		/* Sort the dE's */
+#define HEAP_COMPARE(A,B) (A->de < B->de)
+		HEAPSORT(pval *, stpat, npat);
+#undef HEAP_COMPARE
+
+		/* Create bins and add points */
+		bins[0].min = 0.0;
+		for (nbins = i = 0; i < npat && nbins < maxbins; i++) {
+			double de = stpat[i]->de;
+
+			/* Move on to next bin ? */
+			if (bins[nbins].count >= mincount
+			 && (de - bins[nbins].min) >= mbwidth) {
+				if (i > 0)
+					bins[nbins].max = 0.5 * (de + stpat[i-1]->de);
+				else
+					bins[nbins].max = de;
+				nbins++;
+				bins[nbins].min = bins[nbins-1].max; 
+			} 
+			bins[nbins].count++;
+			bins[nbins].max = de;
+		}
+		if (bins[nbins].count != 0)
+			nbins++;
+
+		/* Compute value */
+		tval = 0.0;
+		for (i = 0; i < nbins; i++) {
+			bins[i].val = bins[i].count/(bins[i].max - bins[i].min);
+			tval += bins[i].val;
+		}
+
+		tval /= 100.0;		/* Make it % */
+		for (i = 0; i < nbins; i++) {
+			bins[i].val /= tval;
+			if (verb) fprintf(verbo,"Bin %d, %f - %f, % 2.4f%%, count %d\n",
+			             i,bins[i].min,bins[i].max,bins[i].val,bins[i].count);
+		}
+
+		/* Plot it */
+		if ((x = (double *)calloc(nbins+1, sizeof(double))) == NULL)
+			error("malloc of histogram x array");
+		if ((y = (double *)calloc(nbins+1, sizeof(double))) == NULL)
+			error("malloc of histogram y array");
+
+		for (i = 0; i < nbins; i++) {
+			x[i] = 0.5 * (bins[i].min + bins[i].max);
+			y[i] = bins[i].val;
+		}
+		x[i] = demax;
+		y[i] = 0.0;
+		do_plot(x, y, NULL, NULL, nbins+1);
+
+		free(y);
+		free(x);
+		free(bins);
+	}
+
+	/* - - - - - - - - - - */
+	/* Create a pruned .ti3 file */
+	if (prune > 0.0) {
+		char *cp, outname[MAXNAMEL+31];
+		cgats *ocg;
+		cgats_set_elem *setel;		/* Array of set value elements */
+
+		strcpy(outname, ti3name);
+		if ((cp = strrchr(outname, '.')) == NULL)
+			cp = outname + strlen(outname);
+		sprintf(cp, "_p%.2f.ti3",prune);
+		
+		if (verb) fprintf(verbo,"Created pruned file '%s'\n",outname);
+
+		/* Create the output files */
+		if ((ocg = new_cgats()) == NULL)
+			error("Failed to create cgats object");
+	
+		/* Duplicate the type of the file */
+		if (icg->t[0].tt == cgats_X) {
+			ocg->add_other(ocg, icg->cgats_type);
+			ocg->add_table(ocg, tt_other, 0);
+		} else if (icg->t[0].tt == tt_other) {
+			ocg->add_other(ocg, icg->others[icg->t[0].oi]);
+			ocg->add_table(ocg, tt_other, 0);
+		} else {
+			ocg->add_table(ocg, icg->t[0].tt, 0);
+		}
+	
+		/* Duplicate all the keywords */
+		for (i = 0; i < icg->t[0].nkwords; i++) {
+			ocg->add_kword(ocg, 0, icg->t[0].ksym[i], icg->t[0].kdata[i], NULL);
+		}
+	
+		/* Duplicate all of the fields */
+		for (i = 0; i < icg->t[0].nfields; i++) {
+			ocg->add_field(ocg, 0, icg->t[0].fsym[i], icg->t[0].ftype[i]);
+		}
+	
+		if ((setel = (cgats_set_elem *)malloc(
+		     sizeof(cgats_set_elem) * icg->t[0].nfields)) == NULL)
+			error("Malloc failed!");
+	
+		/* Copy them approproately */
+		for (i = 0; i < icg->t[0].nsets; i++) {
+
+			if (tpat[i].de <= prune) {
+				icg->get_setarr(icg, 0, i, setel);
+				ocg->add_setarr(ocg, 0, setel);
+			}
+		}
+	
+		if (verb) {
+			double acc;
+
+			acc = (double)ocg->t[0].nsets/(double)icg->t[0].nsets * 100.0;
+			fprintf(verbo,"%.2f%% accepted, %.3f%% rejected\n",acc, 100.0-acc);
+		}
+
+		/* Write out the file */
+		if (ocg->write_name(ocg, outname))
+			error("CGATS file '%s' write error : %s",outname,ocg->err);
+	}
+
+	/* - - - - - - - - - - */
+	/* Display various results */
 	{
 		double merr = 0.0;		/* Max */
 		double aerr = 0.0;		/* Avg */
@@ -853,58 +1097,57 @@ int main(int argc, char *argv[])
 		int inn, outn;			/* Chanells for input and output spaces */
 
 		if (dovrml) {
-			wrl = start_vrml(out_name, doaxes);
-			start_line_set(wrl);
-		}
-
-		/* Open up the file for reading */
-		if ((rd_fp = new_icmFileStd_name(iccname,"r")) == NULL)
-			error("Write: Can't open file '%s'",iccname);
-
-		if ((rd_icco = new_icc()) == NULL)
-			error("Read: Creation of ICC object failed");
-
-		/* Read the header and tag list */
-		if ((rv = rd_icco->read(rd_icco,rd_fp,0)) != 0)
-			error("Read: %d, %s",rv,rd_icco->err);
-
-		/* Get the Fwd table, absolute with Lab override */
-		if ((luo = rd_icco->get_luobj(rd_icco, icmFwd, intent,
-		                              icSigLabData, icmLuOrdNorm)) == NULL) {
-			error("%d, %s",rd_icco->errc, rd_icco->err);
+			wrl = new_vrml(out_name, doaxes, vrml_lab);
+			wrl->start_line_set(wrl, 0);
 		}
 
 		/* Get details of conversion (Arguments may be NULL if info not needed) */
 		luo->spaces(luo, NULL, &inn, NULL, &outn, NULL, NULL, NULL, NULL, NULL);
 
 		for (i = 0; i < npat; i++) {
-			double out[3];
-			double mxd;
+			double de, *out;
 
-			/* Lookup the patch value in the profile */
-			if (luo->lookup(luo, out, tpat[i].p) > 1)
-				error("%d, %s",rd_icco->errc,rd_icco->err);
+			de = tpat[i].de;
+			out = tpat[i].pv;
 
 			if (verb > 1) {
+#ifdef HACK		// Print XYZ
+#pragma message("!!!!!!!!!!!!!!!!!! profile/profcheck.c HACK enabled !!!!!!!!!!!!!!!")
+				double outxyz[3], vxyz[3];
+				icmLab2XYZ(&icmD50, outxyz, out);
+				icmLab2XYZ(&icmD50, vxyz, tpat[i].v);
+
+				printf("[%f] %s%s%s: %s -> %f %f %f should be %f %f %f\n",
+				       de,
+				       tpat[i].sid,
+				       tpat[i].slo[0] != '\000' ? " @ " : "",
+				       tpat[i].slo,
+				       icmPdv(devchan, tpat[i].p),
+				       outxyz[0],outxyz[1],outxyz[2],
+				       vxyz[0],vxyz[1],vxyz[2]);
+#else
 				printf("[%f] %s%s%s: %s -> %f %f %f should be %f %f %f\n",
-				       cie2k ? icmCIE2K(tpat[i].v, out) : 
-				               cie94 ? icmCIE94(tpat[i].v, out) : icmLabDE(tpat[i].v, out),
+				       de,
 				       tpat[i].sid,
 				       tpat[i].slo[0] != '\000' ? " @ " : "",
 				       tpat[i].slo,
 				       icmPdv(devchan, tpat[i].p),
 				       out[0],out[1],out[2],
 				       tpat[i].v[0],tpat[i].v[1],tpat[i].v[2]);
+#endif
 			}
 			if (dovrml) {
-				if (dominl && icmLabDE(tpat[i].v, out) < 0.5) {
+				int ix[2];
+				double p1[3], p2[3];
+
+				/* Add the vertexes */
+				if (dominl && de < 0.5) {		/* Make a minimum length */
 					double cent[3], vec[3], vlen;
-					double p1[3], p2[3];
 
 					/* Compute center */
 					icmAdd3(cent, tpat[i].v, out);
 					icmScale3(cent, cent, 0.5);
-					if ((vlen = icmLabDE(tpat[i].v, out)) < 1e-6) {
+					if ((vlen = de) < 1e-6) {
 						vec[0] = 0.25; vec[1] = 0.0; vec[2] = 0.0;
 					} else {
 						icmSub3(vec, tpat[i].v, out);
@@ -912,36 +1155,36 @@ int main(int argc, char *argv[])
 					}
 					icmSub3(p1, cent, vec);
 					icmAdd3(p2, cent, vec);
-					add_vertex(wrl, p1);
-					add_vertex(wrl, p2);
 				} else {
-					add_vertex(wrl, tpat[i].v);
-					add_vertex(wrl, out);
+					icmCpy3(p1, tpat[i].v);
+					icmCpy3(p2, out);
 				}
-			}
+				ix[0] = wrl->add_vertex(wrl, 0, p1);
+				ix[1] = wrl->add_vertex(wrl, 0, p2);
 
-			/* Check the result */
-			if (cie2k)
-				mxd = icmCIE2K(tpat[i].v, out);
-			else if (cie94)
-				mxd = icmCIE94(tpat[i].v, out);
-			else
-				mxd = icmLabDE(tpat[i].v, out);
+				/* Add the line */
+				if (dodecol) {		/* Lines with color determined by length */
+					double rgb[3];
+					DE2RGB(rgb, icmNorm33(p1, p2));
+					wrl->add_col_line(wrl, 0, ix, rgb);
 
-			aerr += mxd;
-			rerr += mxd * mxd;
+				} else {	/* Natural color */
+					wrl->add_line(wrl, 0, ix);
+				}
+			}
+
+			/* Stats */
+			aerr += de;
+			rerr += de * de;
 
 			nsamps++;
-			if (mxd > merr)
-				merr = mxd;
+			if (de > merr)
+				merr = de;
 
 		}
 		if (dovrml) {
-			if (dodecol)
-				make_de_lines(wrl);
-			else
-				make_lines(wrl, 2);
-			end_vrml(wrl);
+			wrl->make_lines_vc(wrl, 0, 0.0);
+			wrl->del(wrl);
 		}
 		printf("Profile check complete, errors%s: max. = %f, avg. = %f, RMS = %f\n",
             cie2k ? "(CIEDE2000)" : cie94 ? " (CIE94)" : "", merr, aerr/nsamps, sqrt(rerr/nsamps));
@@ -975,13 +1218,13 @@ int main(int argc, char *argv[])
 
 			if (sortbypcs) {
 				/* Sort by pcs delta */
-#define HEAP_COMPARE(A,B) (A.dv < B.dv)
-				HEAPSORT(pval, tpat, npat);
+#define HEAP_COMPARE(A,B) (A->dv < B->dv)
+				HEAPSORT(pval *, stpat, npat);
 #undef HEAP_COMPARE
 			} else {
 				/* Sort by device delta */
-#define HEAP_COMPARE(A,B) (A.dp < B.dp)
-				HEAPSORT(pval, tpat, npat);
+#define HEAP_COMPARE(A,B) (A->dp < B->dp)
+				HEAPSORT(pval *, stpat, npat);
 #undef HEAP_COMPARE
 			}
 
@@ -1000,18 +1243,18 @@ int main(int argc, char *argv[])
 			for (i = 0; i < npat; i++) {
 				if (devspace == icSigCmykData) {
 					printf("%s: %f %f %f %f [%f] -> %f %f %f [%f]\n",
-					       tpat[i].sid,
-					       tpat[i].p[0],tpat[i].p[1],tpat[i].p[2],tpat[i].p[3],
-					       tpat[i].dp,
-					       tpat[i].v[0],tpat[i].v[1],tpat[i].v[2],
-					       tpat[i].dv);
+					       stpat[i]->sid,
+					       stpat[i]->p[0],stpat[i]->p[1],stpat[i]->p[2],stpat[i]->p[3],
+					       stpat[i]->dp,
+					       stpat[i]->v[0],stpat[i]->v[1],stpat[i]->v[2],
+					       stpat[i]->dv);
 				} else {	/* Assume RGB/CMY */
 					printf("%s: %f %f %f [%f] -> %f %f %f [%f]\n",
-					       tpat[i].sid,
-					       tpat[i].p[0],tpat[i].p[1],tpat[i].p[2],
-					       tpat[i].dp,
-					       tpat[i].v[0],tpat[i].v[1],tpat[i].v[2],
-					       tpat[i].dv);
+					       stpat[i]->sid,
+					       stpat[i]->p[0],stpat[i]->p[1],stpat[i]->p[2],
+					       stpat[i]->dp,
+					       stpat[i]->v[0],stpat[i]->v[1],stpat[i]->v[2],
+					       stpat[i]->dv);
 				}
 			}
 		}
@@ -1024,292 +1267,13 @@ int main(int argc, char *argv[])
 		rd_fp->del(rd_fp);
 	}
 
+	icg->del(icg);		/* Clean up */
+
 	return 0;
 }
 
 
 /* ------------------------------------------------ */
-/* Some simple functions to do basix VRML work */
-/* !!! Should change to plot/vrml lib !!! */
-
-#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);
-static void DE2RGB(double *out, double in);
-
-FILE *start_vrml(char *name, int doaxes) {
-	FILE *wrl;
-
-	/* Define the axis boxes */
-	struct {
-		double x, y, z;			/* Box center */
-		double wx, wy, wz;		/* Box size */
-		double r, g, b;			/* Box color */
-	} 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 */
-	};
-
-	/* Define the labels */
-	struct {
-		double x, y, z;
-		double size;
-		char *string;
-		double r, g, b;
-	} labels[6] = {
-		{ -2, 2, -GAMUT_LCENT + 100 + 10, 10, "+L*",  .7, .7, .7 },	/* Top of L axis */
-		{ -2, 2, -GAMUT_LCENT - 10,      10, "0",    .7, .7, .7 },	/* Bottom of L axis */
-		{ 100 + 5, -3,  0-GAMUT_LCENT,  10, "+a*",  1,  0,  0 },	/* +a (red) axis */
-		{ -5, -100 - 10, 0-GAMUT_LCENT,  10, "-b*",  0,  0,  1 },	/* -b (blue) axis */
-		{ -100 - 15, -3, 0-GAMUT_LCENT,  10, "-a*",  0,  0,  1 },	/* -a (green) axis */
-		{ -5,  100 + 5, 0-GAMUT_LCENT,  10, "+b*",  1,  1,  0 },	/* +b (yellow) axis */
-	};
-
-	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) {
-		int n;
-		fprintf(wrl,"    # Lab axes as boxes:\n");
-		for (n = 0; n < 5; n++) {
-			fprintf(wrl,"    Transform { translation %f %f %f\n", axes[n].x, axes[n].y, axes[n].z);
-			fprintf(wrl,"      children [\n");
-			fprintf(wrl,"        Shape{\n");
-			fprintf(wrl,"          geometry Box { size %f %f %f }\n",
-			                       axes[n].wx, axes[n].wy, axes[n].wz);
-			fprintf(wrl,"          appearance Appearance { material Material ");
-			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[n].r, axes[n].g, axes[n].b);
-			fprintf(wrl,"        }\n");
-			fprintf(wrl,"      ]\n");
-			fprintf(wrl,"    }\n");
-		}
-		fprintf(wrl,"    # Axes identification:\n");
-		for (n = 0; n < 6; n++) {
-			fprintf(wrl,"    Transform { translation %f %f %f\n", labels[n].x, labels[n].y, labels[n].z);
-			fprintf(wrl,"      children [\n");
-			fprintf(wrl,"        Shape{\n");
-			fprintf(wrl,"          geometry Text { string [\"%s\"]\n",labels[n].string);
-			fprintf(wrl,"            fontStyle FontStyle { family \"SANS\" style \"BOLD\" size %f }\n",
-			                                  labels[n].size);
-			fprintf(wrl,"                        }\n");
-			fprintf(wrl,"          appearance Appearance { material Material ");
-			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", labels[n].r, labels[n].g, labels[n].b);
-			fprintf(wrl,"        }\n");
-			fprintf(wrl,"      ]\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");
-}
-
-/* Assume 2 ppset, and make line color prop to length */
-void make_de_lines(FILE *wrl) {
-	int i, j;
-
-	fprintf(wrl,"      ]\n");
-	fprintf(wrl,"    }\n");
-	fprintf(wrl,"  coordIndex [\n");
-
-	for (i = 0; i < npoints;) {
-		for (j = 0; j < 2; 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], ss;
-		for (ss = 0.0, j = 0; j < 3; j++) {
-			double tt = (pary[i & ~1].pp[j] - pary[i | 1].pp[j]);
-			ss += tt * tt;
-		}
-		ss = sqrt(ss);
-		DE2RGB(rgb, ss);
-		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;
-}
 
 /* Convert a delta E value into a signal color: */
 static void
diff --git a/profile/profin.c b/profile/profin.c
index c703c42..98cad65 100644
--- a/profile/profin.c
+++ b/profile/profin.c
@@ -678,15 +678,15 @@ make_input_icc(
 		if ((ri = icg->find_field(icg, 0, "RGB_R")) < 0)
 			error ("Input file doesn't contain field RGB_R");
 		if (icg->t[0].ftype[ri] != r_t)
-			error ("Field CMYK_C is wrong type - corrupted file ?");
+			error ("Field RGB_R is wrong type - corrupted file ?");
 		if ((gi = icg->find_field(icg, 0, "RGB_G")) < 0)
 			error ("Input file doesn't contain field RGB_G");
 		if (icg->t[0].ftype[gi] != r_t)
-			error ("Field CMYK_M is wrong type - corrupted file ?");
+			error ("Field RGB_G is wrong type - corrupted file ?");
 		if ((bi = icg->find_field(icg, 0, "RGB_B")) < 0)
 			error ("Input file doesn't contain field RGB_B");
 		if (icg->t[0].ftype[bi] != r_t)
-			error ("Field CMYK_Y is wrong type - corrupted file ?");
+			error ("Field RGB_B is wrong type - corrupted file ?");
 
 		if (spec == 0) {        /* Using instrument tristimulous value */
 
@@ -838,7 +838,9 @@ make_input_icc(
 		               wr_xicc, icmFwd, icmDefaultIntent,
 		               icmLuOrdNorm,
 		               flags, 		/* Flags */
-		               npat, npat, tpat, NULL, 0.0, wpscale, smooth, avgdev, 1.0,
+		               npat, npat, tpat, NULL, 0.0, wpscale,
+//		               NULL,		/* bpo */
+			           smooth, avgdev, 1.0,
 		               NULL, NULL, NULL, iquality)) == NULL)
 			error("%d, %s",wr_xicc->errc, wr_xicc->err);
 
@@ -890,7 +892,7 @@ make_input_icc(
 			/* Create gamma/matrix model to extrapolate with. */
 			/* (Use ofset & gain, gamma curve as 0th order with 1 harmonic, */
 			/* and smooth it.) */
-			if ((mm = new_MatrixModel(verb, nmpat, mpat, wantLab,
+			if ((mm = new_MatrixModel(wr_icco, verb, nmpat, mpat, wantLab,
 				      /* quality */ -1, /* isLinear */ ptype == prof_matonly,
 				      /* isGamma */ 0, /* isShTRC */ 0,
 				      /* shape0gam */ 1, /* clipbw */ 0, /* clipprims */ 0,
@@ -1123,7 +1125,9 @@ make_input_icc(
 			               ICX_2PASSSMTH |
 #endif
 		               flags, 		/* Flags */
-		               npat + nxpat, npat, tpat, NULL, 0.0, wpscale, smooth, avgdev, 1.0,
+		               npat + nxpat, npat, tpat, NULL, 0.0, wpscale,
+//			           NULL,	/* bpo */
+			           smooth, avgdev, 1.0,
 			           NULL, NULL, NULL, iquality)) == NULL)
 			error ("%d, %s",wr_xicc->errc, wr_xicc->err);
 
diff --git a/profile/profout.c b/profile/profout.c
index 8c6fa0d..790ffa4 100644
--- a/profile/profout.c
+++ b/profile/profout.c
@@ -66,6 +66,7 @@
 #undef IMP_MONO					/* [Undef] Turn on development code */
 
 #define EMPH_DISP_BLACKPOINT	/* [Define] Increase weight near diplay black point */
+#define IGNORE_DISP_ZEROS	    /* [Define] Ignore points with zero value if not at dev. zero */
 #define NO_B2A_PCS_CURVES		/* [Define] PCS curves seem to make B2A less accurate. Why ? */
 #define USE_CAM_CLIP_OPT		/* [Define] Clip out of gamut in CAM space rather than PCS */
 #define USE_LEASTSQUARES_APROX	/* [Define] Use least squares fitting approximation in B2A */
@@ -126,7 +127,7 @@
 	looking up the profile in absolute mode.
 
 	For a Matrix profile, in the case of the white point this is
-	because we're not using the ICC 16 bit quantized value to
+	because we may not be using the ICC 16 bit quantized value to
 	create the relative transform matrix, and in the case of
 	the black point, it can never be a perfect match because the black
 	point returned by a profile lookup will be the quantized black
@@ -644,6 +645,7 @@ make_output_icc(
 	int verify,				/* nz to print verification */
 	int clipprims,			/* Clip white, black and primaries */
 	double wpscale,			/* >= 0.0 for media white point scale factor */
+//	double *bpo,			/* != NULL for XYZ black point override */
 	icxInk *oink,			/* Ink limit/black generation setup (NULL if n/a) */
 	char *in_name,			/* input .ti3 file name */
 	char *file_name,		/* output icc name */
@@ -671,8 +673,8 @@ make_output_icc(
 	profxinf *xpi			/* Optional Profile creation extra data */
 ) {
 	int isdisp;				/* nz if this is a display device, 0 if output */
-	double dispLuminance = 0.0;	/* Display luminance. 0 if not known */
-	int isdnormed = 0;		/* Has display data been normalised to 100 ? */
+	double dispLuminance = 0.0;	/* Display luminance. 0.0 if not known */
+	int isdnormed = 0;		/* Has display data been normalised to white Y = 100 ? */
 	int allintents;			/* nz if all intents should possibly be created */
 	icmFile *wr_fp;
 	icc *wr_icco;
@@ -745,6 +747,8 @@ make_output_icc(
 		}
 
 		/* See if the CIE data has been normalised to Y = 100 */
+		/* If so, it's assumed to be by LUMINANCE_XYZ_CDM2 */
+		/* By default assume not. */
 		if ((ti = icg->find_kword(icg, 0, "NORMALIZED_TO_Y_100")) < 0
 		 || strcmp(icg->t[0].kdata[ti],"NO") == 0) {
 			isdnormed = 0;
@@ -1020,15 +1024,7 @@ make_output_icc(
 		if (xpi != NULL && xpi->creator != 0L)
 			wh->creator = xpi->creator;
 
-#ifdef NT
-		wh->platform = icSigMicrosoft;
-#endif
-#ifdef __APPLE__
-		wh->platform = icSigMacintosh;
-#endif
-#if defined(UNIX) && !defined(__APPLE__)
-		wh->platform = icmSig_nix;
-#endif
+		/* Default platform is OK */
 
 		if (xpi != NULL && xpi->transparency)
 			wh->attributes.l |= icTransparency;
@@ -1656,15 +1652,34 @@ make_output_icc(
 				tpat[i].v[0] = *((double *)icg->t[0].fdata[i][Xi]);
 				tpat[i].v[1] = *((double *)icg->t[0].fdata[i][Yi]);
 				tpat[i].v[2] = *((double *)icg->t[0].fdata[i][Zi]);
-				if (!isLab && (!isdisp || isdnormed != 0)) {
-					tpat[i].v[0] /= 100.0;		/* Normalise XYZ to range 0.0 - 1.0 */
-					tpat[i].v[1] /= 100.0;
-					tpat[i].v[2] /= 100.0;
-				}
-				if (!isLab && wantLab) { /* Convert test patch result XYZ to PCS (D50 Lab) */
-					icmXYZ2Lab(&icmD50, tpat[i].v, tpat[i].v);
-				} else if (isLab && !wantLab) {
-					icmLab2XYZ(&icmD50, tpat[i].v, tpat[i].v);
+				/* For display, convert to measurement XYZ and re-normalise later */
+				if (isdisp) {
+					if (isLab) {
+						icmLab2XYZ(&icmD50, tpat[i].v, tpat[i].v);
+						isLab = 0;
+					} else if (isdnormed) {
+						tpat[i].v[0] /= 100.0;		/* Normalise XYZ to range 0.0 - 1.0 */
+						tpat[i].v[1] /= 100.0;
+						tpat[i].v[2] /= 100.0;
+					}
+					if (isdnormed && dispLuminance > 0.0) {
+						tpat[i].v[0] *= dispLuminance;
+						tpat[i].v[1] *= dispLuminance;
+						tpat[i].v[2] *= dispLuminance;
+					}	/* else Hmm. */
+
+				/* Convert to normalised 0.0 - 1.0 range in target PCS */
+				} else {
+					if (!isLab) {
+						tpat[i].v[0] /= 100.0;		/* Normalise XYZ to range 0.0 - 1.0 */
+						tpat[i].v[1] /= 100.0;
+						tpat[i].v[2] /= 100.0;
+					}
+					if (!isLab && wantLab) { /* Convert test patch result XYZ to PCS (D50 Lab) */
+						icmXYZ2Lab(&icmD50, tpat[i].v, tpat[i].v);
+					} else if (isLab && !wantLab) {
+						icmLab2XYZ(&icmD50, tpat[i].v, tpat[i].v);
+					}
 				}
 			}
 
@@ -1684,10 +1699,14 @@ make_output_icc(
 			if ((ii = icg->find_kword(icg, 0, "SPECTRAL_END_NM")) < 0)
 				error ("Input file doesn't contain keyword SPECTRAL_END_NM");
 			sp.spec_wl_long = atof(icg->t[0].kdata[ii]);
-			if (!isdisp || isdnormed != 0)
-				sp.norm = 100.0;
-			else
+			if (isdisp) {				/* convert to measurement values - re-norm later */
 				sp.norm = 1.0;
+				if (isdnormed)
+					sp.norm *= 100.0;
+				if (isdnormed && dispLuminance > 0.0)
+					sp.norm /= dispLuminance;
+			} else
+				sp.norm = 100.0;		/* Convert to 0.0 - 1.0 ref/trans range */
 
 			/* Find the fields for spectral values */
 			for (j = 0; j < sp.spec_n; j++) {
@@ -1820,7 +1839,6 @@ make_output_icc(
 			}
 
 			for (i = 0; i < npat; i++) {
-
 				tpat[i].w = 1.0;
 				tpat[i].p[0] = *((double *)icg->t[0].fdata[i][ci]) / 100.0;
 				tpat[i].p[1] = *((double *)icg->t[0].fdata[i][mi]) / 100.0;
@@ -1841,59 +1859,45 @@ make_output_icc(
 
 				/* Convert it to CIE space */
 				sp2cie->convert(sp2cie, tpat[i].v, &sp);
-
 			}
 
 			sp2cie->del(sp2cie);		/* Done with this */
-
 		}
 
-		isLab = wantLab;		/* We now have what we want */
-
-		/* Normalize display values to Y = 1.0 if needed */
-		/* (re-norm spec derived, since observer may be different) */
-		if (isdisp && (isdnormed == 0 || spec != 0)) {
+		/* Normalize display values to Y = 1.0 */
+		if (isdisp) {
 			double scale = -1e6;
 
-			if (wantLab) {
-				double bxyz[3];
+			if (isLab)		/* assert */
+				error("Internal - display values must be XYZ for normalisation");
 
-				/* Locate max Y */
-				for (i = 0; i < npat; i++) {
-					icmLab2XYZ(&icmD50, bxyz,  tpat[i].v);
-					if (bxyz[1] > scale)
-						scale = bxyz[1];
-				}
-				
-				scale = 1.0/scale;
-
-				/* Scale max Y to 1.0 */
-				for (i = 0; i < npat; i++) {
-					icmLab2XYZ(&icmD50, tpat[i].v, tpat[i].v);
-					tpat[i].v[0] *= scale;
-					tpat[i].v[1] *= scale;
-					tpat[i].v[2] *= scale;
-					icmXYZ2Lab(&icmD50, tpat[i].v, tpat[i].v);
-				}
-			} else {
+			/* Locate max Y */
+			for (i = 0; i < npat; i++) {
+				if (tpat[i].v[1] > scale)
+					scale = tpat[i].v[1];
+			}
+			
+			scale = 1.0/scale;
 
-				/* Locate max Y */
-				for (i = 0; i < npat; i++) {
-					if (tpat[i].v[1] > scale)
-						scale = tpat[i].v[1];
-				}
-				
-				scale = 1.0/scale;
+			for (i = 0; i < npat; i++) {
+				tpat[i].v[0] *= scale;
+				tpat[i].v[1] *= scale;
+				tpat[i].v[2] *= scale;
 
-				for (i = 0; i < npat; i++) {
-					tpat[i].v[0] *= scale;
-					tpat[i].v[1] *= scale;
-					tpat[i].v[2] *= scale;
-				}
+				if (wantLab)
+					icmXYZ2Lab(&icmD50, tpat[i].v, tpat[i].v);
 			}
+
+			/* Change Black Point Override XYZ from reading to normalised */
+//			if (bpo != NULL) {
+//				bpo[0] *= scale;
+//				bpo[1] *= scale;
+//				bpo[2] *= scale;
+//			}
 		}
-	}	/* End of reading in CGATs file */
+		isLab = wantLab;		/* We now have what we want */
 
+	}	/* End of reading in CGATs file */
 
 #ifdef EMPH_DISP_BLACKPOINT
 	/* Add extra weighting to sample points near black for additive display. */
@@ -1903,20 +1907,20 @@ make_output_icc(
 	/* of the possible "scaled" viewing mode of additive display */
 	/* usage ? What about print and scan ?? */
 	if (isdisp && (imask == ICX_W || imask == ICX_RGB)) {
-		double minL = 1e6;;
+		double minL = 1e6;
 
 		/* Locate the lowest L* value */
 		for (i = 0; i < npat; i++) {
-			if (wantLab) {
-				if (tpat[i].v[0] < minL)
-					minL = tpat[i].v[0];
-			} else {
-				double lab[2];
+			double lab[3];
+
+			if (wantLab)
+				lab[0] = tpat[i].v[0];
+			else
 				icmXYZ2Lab(&icmD50, lab, tpat[i].v);
-				if (lab[0] < minL)
-					minL = lab[0];
-			}
+			if (lab[0] < minL)
+				minL = lab[0];
 		}
+//printf("~1 final minL = %f\n",minL);
 
 		/* Compute weighting factor */
 		/* (Hard to guess what numbers to put in here.. */
@@ -1934,11 +1938,52 @@ make_output_icc(
 				continue;
 			L = 1.0 + 19.0 * pow(1.0 - L, 3.0); 
 			tpat[i].w *= L;
-//printf("~1 pat %d %f %f %f weight %f\n", i,tpat[i].p[0], tpat[i].p[1], tpat[i].p[2], tpat[i].w);
+//printf("~1 pat %d: %f %f %f weight %f\n", i,tpat[i].p[0], tpat[i].p[1], tpat[i].p[2], tpat[i].w);
 		}
 	}
 #endif /* EMPH_DISP_BLACKPOINT */
 
+#ifdef IGNORE_DISP_ZEROS
+	/* If a display has a very good black, and the instrument is not sensitive */
+	/* enough to properly measur the near black values and returns 0.0, */
+	/* then the resulting profile will tend to incorrectly boost the */
+	/* dark shadows. A heuristic to counteract this problem is to */
+	/* ignore any readings that have any value <= 0.0 except */
+	/* those for device black. */ 
+	
+	if (isdisp && (imask == ICX_W || imask == ICX_RGB)) {
+		int noomit = 0;
+
+		for (i = 0; i < npat; i++) {
+			double xyz[3], L;
+			if (wantLab) 
+				icmLab2XYZ(&icmD50, xyz, tpat[i].v);
+			else
+				icmCpy3(xyz, tpat[i].v);
+
+			/* Don't ignore device zero point */
+			if (imask == ICX_W) {
+				if (tpat[i].p[0] <= 0.0)
+					continue;
+			} else {
+				if (tpat[i].p[0] <= 0.0
+				 && tpat[i].p[1] <= 0.0
+				 && tpat[i].p[2] <= 0.0)
+					continue;
+			}
+//printf("~1 pat %d: XYZ %f %f %f\n", i, xyz[0], xyz[1], xyz[2]);
+			/* Ignore any XYZ that is zero */
+			if (xyz[0] <= 0.0 || xyz[1] <= 0.0 || xyz[2] <= 0.0) {
+				tpat[i].w = 0.0;
+				noomit++;
+//printf("~1  ignored\n");
+			}
+		}
+		if (verb)
+			fprintf(verbo,"Omitted %d zero measurements\n",noomit);
+	}
+#endif /* IGNORE_DISP_ZEROS */
+
 	if (isLut) {
 		xicc *wr_xicc;			/* extention object */
 		icxLuBase *AtoB;		/* AtoB ixcLu */
@@ -1981,7 +2026,9 @@ make_output_icc(
 			               ICX_2PASSSMTH |
 #endif
 			               flags,
-			               npat, npat, tpat, NULL, dispLuminance, wpscale, smooth, avgdev, demph, 
+			               npat, npat, tpat, NULL, dispLuminance, wpscale,
+//				           bpo,
+				           smooth, avgdev, demph, 
 			               NULL, oink, cal, iquality)) == NULL)
 				error("%d, %s",wr_xicc->errc, wr_xicc->err);
 
@@ -2187,6 +2234,9 @@ make_output_icc(
 						if (absname[i] == absname[j]) {
 							cx.abs_intent[j] = cx.abs_intent[i];
 							cx.abs_luo[j] = cx.abs_luo[i];
+							abs_xicc[j] = abs_xicc[i];
+							abs_icc[j] = abs_icc[i];
+							abs_fp[j] = abs_fp[i];
 							if (verb)
 								printf("Applying %s abstract profile '%s' to %s table\n",
 								i == 0 ? "first" : i == 1 ? "second" : "third",
@@ -2490,7 +2540,7 @@ make_output_icc(
 
 				optcomb tcomb = oc_imo;	/* Create all by default */
 
-				if ((xf = new_xfit()) == NULL) {
+				if ((xf = new_xfit(icco)) == NULL) {
 					error("profout: Creation of xfit object failed");
 				}
 					
@@ -2727,15 +2777,17 @@ make_output_icc(
 			/* Free up abstract transform */
 			for (i = 0; i < cx.ntables; i++) {
 				if (cx.abs_luo[i] != NULL) {
-					for (j = cx.ntables-1; j >= i; j--) {	/* Free all duplicates */
-						if (cx.abs_luo[j] == cx.abs_luo[i]) {
-							cx.abs_luo[j]->del(cx.abs_luo[j]);
-							abs_xicc[j]->del(abs_xicc[j]);
-							abs_icc[j]->del(abs_icc[j]);
-							abs_fp[j]->del(abs_fp[j]);
+					/* Free this and all associated resources */
+					cx.abs_luo[i]->del(cx.abs_luo[i]);
+					abs_xicc[i]->del(abs_xicc[i]);
+					abs_icc[i]->del(abs_icc[i]);
+					abs_fp[i]->del(abs_fp[i]);
+					/* Mark all duplicates as being free'd too */
+					for (j = i+1; j < cx.ntables; j++) {
+						if (cx.abs_luo[j] == cx.abs_luo[i])
 							cx.abs_luo[j] = NULL;
-						}
 					}
+					cx.abs_luo[i] = NULL;
 				}
 			}
 
@@ -2899,7 +2951,9 @@ make_output_icc(
 		               wr_xicc, icmFwd, isdisp ? icmDefaultIntent : icRelativeColorimetric,
 		               icmLuOrdRev,
 			           flags,		 		/* Compute white & black */
-		               npat, npat, tpat, NULL, dispLuminance, wpscale, smooth, avgdev, demph,
+		               npat, npat, tpat, NULL, dispLuminance, wpscale,
+//			           bpo,
+			           smooth, avgdev, demph,
 		               NULL, oink, cal, iquality)) == NULL)
 			error("%d, %s",wr_xicc->errc, wr_xicc->err);
 
diff --git a/profile/splitti3.c b/profile/splitti3.c
index 0e170c1..a48c6fa 100644
--- a/profile/splitti3.c
+++ b/profile/splitti3.c
@@ -50,7 +50,7 @@ void
 usage(void) {
 	fprintf(stderr,"Split a .ti3 into two, Version %s\n",ARGYLL_VERSION_STR);
 	fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
-	fprintf(stderr,"usage: splitcgats [-options] input.ti3 output1.ti3 output2.ti3\n");
+	fprintf(stderr,"usage: splitti3 [-options] input.ti3 output1.ti3 output2.ti3\n");
 	fprintf(stderr," -v              Verbose - print each patch value\n");
 	fprintf(stderr," -n no           Put no sets in first file, and balance in second file.\n");
 	fprintf(stderr," -p percent      Put percent%% sets in first file, and balance in second file. (def. 50%%)\n");
diff --git a/profile/txt2ti3.c b/profile/txt2ti3.c
index 62807f3..c1d9ed0 100644
--- a/profile/txt2ti3.c
+++ b/profile/txt2ti3.c
@@ -43,6 +43,7 @@
 #include "xspect.h"
 #include "insttypes.h"
 #include "numlib.h"
+#include "ui.h"
 
 void
 usage(char *mes) {
@@ -200,7 +201,7 @@ int main(int argc, char *argv[])
 	cmy->add_other(cmy, "ECI2002");			/* Gretag/Logo Target file */
 	cmy->add_other(cmy, ""); 				/* Wildcard */
 	if (cmy->read_name(cmy, devname))
-		error ("Read: Can't read dev file '%s'. Unknown format or corrupted file ? (%s)",devname,cmy->err);
+		error ("Read: Can't read dev file '%s'. Unknown format, missing or corrupted file ? (%s)",devname,cmy->err);
 	if (cmy->ntables != 1)
 		warning("Input file '%s' doesn't contain exactly one table",devname);
 
@@ -291,7 +292,7 @@ int main(int argc, char *argv[])
 	ncie->add_other(ncie, "ECI2002"); 		/* Gretag/Logo Target file */
 	ncie->add_other(ncie, ""); 				/* Wildcard */
 	if (ncie->read_name(ncie, ciename))
-		error ("Read: Can't read cie file '%s'. Unknown format or corrupted file ?",ciename);
+		error ("Read: Can't read cie file '%s'. Unknown format, missing or corrupted file ?",ciename);
 	if (ncie->ntables != 1)
 		warning("Input file '%s' doesn't contain exactly one table",ciename);
 
@@ -365,7 +366,7 @@ int main(int argc, char *argv[])
 		spec->add_other(spec, "ECI2002"); 		/* Gretag/Logo Target file */
 		spec->add_other(spec, ""); 				/* Wildcard */
 		if (spec->read_name(spec, specname))
-			error ("Read: Can't read spec file '%s'. Unknown format or corrupted file ?",specname);
+			error ("Read: Can't read spec file '%s'. Unknown format, missing or corrupted file ?",specname);
 		if (spec->ntables != 1)
 			warning("Input file '%s' doesn't contain exactly one table",specname);