Merge branch 'feature/upstream' into develop

1 files changed, 797 insertions, 166 deletions

diff --git a/rspl/rev.c b/rspl/rev.c
index f2f717b..017943b 100644..100755
--- a/rspl/rev.c
+++ b/rspl/rev.c
@@ -72,7 +72,6 @@
 	a particular aux target that lies within that segment.
 	(1150 near black, k ~= 0.4).
 
-
  */
 
 #include <stdio.h>
@@ -109,9 +108,13 @@
 //#undef DMALLOC_GLOBALS
 
 #define DOSORT			/* [def] Cell sort for better speed */
+#undef EN_UNTWIST		/* [und] Force attempt to try and untwist gamut surface */
+						/* - Seems to improve some, make some worse ?? (i.e Bonet) */
+						/* By default this is controlled using ARGYLL_UNTWIST_GAMUT_SURFACE */
+						/* environment variable. */
 
 #undef REVTABLESTATS	/* [und] Reverse table stats */
-#undef REVVRML			/* [und] Reverse table plots */
+#undef REVVRML			/* [und] Reverse table 3D plots */
 
 #undef DEBUG1			/* [und] Higher level code */
 #undef DEBUG2			/* [und] Lower level code */
@@ -172,16 +175,16 @@ int thissz, lastsz = -1;
 #endif
 
 #ifdef REVTABLESTATS
-#pragma message("!!!!!!!!! REVTABLESTATS set in rspl/rev.c !!!!!!!!!!!")
+# pragma message("!!!!!!!!! REVTABLESTATS set in rspl/rev.c !!!!!!!!!!!")
 #endif
 
 #ifdef REVVRML
-#pragma message("!!!!!!!!! REVVRML set in rspl/rev.c !!!!!!!!!!!")
-#include "vrml.h"
+# pragma message("!!!!!!!!! REVVRML set in rspl/rev.c !!!!!!!!!!!")
+# include "vrml.h"
 #endif
 
 #ifdef CHECK_NNLU
-#pragma message("!!!!!!!!! CHECK_NNLU set in rspl/rspl.h !!!!!!!!!!!")
+# pragma message("!!!!!!!!! CHECK_NNLU set in rspl/rspl.h !!!!!!!!!!!")
 #endif
 
 /* Do an arbitrary printf */
@@ -356,7 +359,7 @@ static void rev_test_vram(size_t size) {
 static void *rev_malloc(rspl *s, size_t size) {
 	void *rv;
 
-	if ((size + 1 * 1024 * 1204) > g_test_ram)
+	if ((size + 1 * 1024 * 1024) > g_test_ram)
 		rev_test_vram(size);
 	if ((rv = malloc(size)) == NULL) {
 		rev_reduce_cache(size);
@@ -371,7 +374,7 @@ static void *rev_malloc(rspl *s, size_t size) {
 static void *rev_calloc(rspl *s, size_t num, size_t size) {
 	void *rv;
 
-	if (((num * size) + 1 * 1024 * 1204) > g_test_ram)
+	if (((num * size) + 1 * 1024 * 1024) > g_test_ram)
 		rev_test_vram(size);
 	if ((rv = calloc(num, size)) == NULL) {
 		rev_reduce_cache(num * size);
@@ -386,7 +389,7 @@ static void *rev_calloc(rspl *s, size_t num, size_t size) {
 static void *rev_realloc(rspl *s, void *ptr, size_t size) {
 	void *rv;
 
-	if ((size + 1 * 1024 * 1204) > g_test_ram)
+	if ((size + 1 * 1024 * 1024) > g_test_ram)
 		rev_test_vram(size);
 	if ((rv = realloc(ptr, size)) == NULL) {
 		rev_reduce_cache(size);		/* approximation */
@@ -1323,7 +1326,7 @@ set_search_limit(
 					/* should evaluate in[0..di-1], and return number that is not to exceed */
 					/* limitv. NULL if not used */
 	void *lcntx,	/* Context passed to limit() */
-	double limitv	/* Value that limit() is not to exceed */
+	double limitv	/* Value that limitf() is not to exceed */
 ) {
 	schbase *b = NULL;		/* Pointer to search base information structure */
 
@@ -4025,10 +4028,11 @@ double *err		/* Output error (weighted) distance at solution point */
 				}
 
 #ifdef DEBUG
-				if (wsrv2)
+				if (wsrv2) {
 					DBG(("Got second ink limit triangle in tetrahedron\n"));
-				else if (wsrv)
+				} else if (wsrv) {
 					DBG(("Got first ink limit triangle in tetrahedron\n"));
+				}
 #endif
 				*err = dist;
 				return wsrv;
@@ -6820,7 +6824,7 @@ static void free_surfhash(rspl *s, int del) {
 
 /* vertex status */
 typedef enum {
-	vtx_norm = 0,		/* Normal vertex in primary bxcell */
+	vtx_norm = 0,		/* Normal vertex in primary bxcell - initial value */
 	vtx_sha  = 1,		/* Vertex has been shadowed */
 	vtx_del  = 2,		/* Vertex has been deleted because it's shadowed */
 	vtx_oil  = 3		/* Vertex is over ink limit */
@@ -6828,8 +6832,9 @@ typedef enum {
 
 struct _vtxrec {
 	int ix;						/* fwd index of vertex */
-	int cix;					/* Cell index for this vertex */
-	double v[MXRO];				/* Output value of vertex */
+	int cix;					/* Fwd cell vertex is in index */
+	double vv[MXRO];			/* Output value of vertex */
+	double vl[MXRO];			/* Log compressed output value of vertex */
 	double dist;				/* Distance from center point squared */
 	int tcount;					/* Touch count for converting to fwd cells */
 	int acount;					/* Actual count for converting to fwd cells */
@@ -6866,7 +6871,7 @@ struct _vtxcache {
 }; typedef struct _vtxcache vtxcache;
 
 
-/* Create the vertex list & hash */
+/* Create the (empty) vertex list & hash */
 static void create_vtxrec_list(rspl *s, vtxcache *vc) {
 	vc->hash_size = primes[3];		/* 3373 */
 	if ((vc->hash = (vtxrec **) rev_calloc(s, vc->hash_size, sizeof(vtxrec *))) == NULL)
@@ -6949,16 +6954,67 @@ static vtxrec *get_vtxrec(vtxcache *vc, int ix) {
 
 	hash = ix % vc->hash_size;
 
-	for (vx = vc->hash[hash]; vx != NULL; vx = vx->hlink) {
+	for (vx = vc->hash[hash]; vx != NULL; vx = vx->hlink) { 
 		if (vx->ix == ix)
 			return vx;
 	}
 	return NULL;
 }
 
+/* Log compress an output value wrt to center point */
+static void logcomp(
+	rspl *s,
+	double *out,
+	double *in,
+	double *cent
+) {
+	int f, fdi = s->fdi;
+	double len;
+
+	if (s->rev.surflin_en) {
+#ifdef NEVER
+		/* (This doesn't seem to improve things) */
+		/* Calculate vector length */
+		for (len = 0.0, f = 0; f < fdi; f++) {
+			double tt= in[f] - cent[f]; 
+			len += tt * tt;
+		}
+		len = sqrt(len);
+
+		/* change length to log length */
+		if (len > DBL_EPSILON) {
+ 
+			len = 20.0 * pow(len, 0.25)/len;	/* Ratio */ 
+
+			for (f = 0; f < fdi; f++) {
+				double tt = in[f] - cent[f]; 
+				out[f] = len * tt + cent[f];
+			}
+		}
+#else
+		if (s->rev.surflin != NULL) {
+			co p;
+
+			for (f = 0; f < fdi; f++)
+				p.p[f] = in[f];
+			s->rev.surflin->interp(s->rev.surflin, &p);
+			for (f = 0; f < fdi; f++)
+				out[f] = p.v[f] - s->rev.linoff[f];
+		} else {
+			for (f = 0; f < fdi; f++)
+				out[f] = in[f];
+		}
+#endif
+	} else {
+		for (f = 0; f < fdi; f++)
+			out[f] = in[f];
+	}
+}
+
 /* Create a new vtxrec or return the current one. */
 /* Allocates it, adds it to cache. */
 /* DOESN"T add it to vtxlist. */
+/* If new, sets status = vtx_norm */
 static vtxrec *new_vtxrec(
 	rspl *s,
 	vtxcache *vc,
@@ -7005,12 +7061,15 @@ static vtxrec *new_vtxrec(
 
 	/* Get the output value */
 	for (f = 0; f < fdi; f++)
-		vx->v[f] = gp[f];
+		vx->vv[f] = gp[f];
+
+	/* Set vl[] */
+	logcomp(s, vx->vl, vx->vv, s->rev.ocent);
 
 	/* Compute distance to overall center point squared */
 	vx->dist = 0.0;
 	for (f = 0; f < fdi; f++) {
-		double tt = gp[f] - s->rev.ocent[f]; 
+		double tt = vx->vl[f] - s->rev.ocent[f]; 
 		vx->dist += tt * tt;
 	}
 
@@ -7020,7 +7079,7 @@ static vtxrec *new_vtxrec(
 		int mi;
 		double gw = s->rev.gw[f];
 		double gl = s->rev.gl[f];
-		t = (vx->v[f] - gl)/gw;
+		t = (vx->vv[f] - gl)/gw;
 		mi = (int)floor(t);			/* Grid coordinate */
 		if (mi < 0)					/* Limit to valid cube base index range */
 			mi = 0;
@@ -7442,7 +7501,7 @@ static void extend_bxcell_shadow_group(
 
 		/* vertex length from gamut center */
 		for (vlen= 0.0, f = 0; f < fdi; f++) {
-			double tt = vx->v[f] - s->rev.ocent[f];
+			double tt = vx->vl[f] - s->rev.ocent[f];
 			vlen += tt * tt;
 		}
 		vlen = sqrt(vlen);
@@ -7453,7 +7512,7 @@ static void extend_bxcell_shadow_group(
 			scale = 1.0;
 
 		for (f = 0; f < fdi; f++)
-			sv[f] = (scale * (vx->v[f] - s->rev.ocent[f])) + s->rev.ocent[f];
+			sv[f] = (scale * (vx->vl[f] - s->rev.ocent[f])) + s->rev.ocent[f];
 
 		/* Distance from scaled vertex to group center */
 		for (w = 0.0, f = 0; f < fdi; f++) {
@@ -8389,13 +8448,13 @@ static void plot_fxcell_surface(rspl *s, int dofclabels, int dobxcells, int dowa
 	is inside or on the surface of the gamut. A simple and definitive
 	topological rule hasn't been forthcoming, so a simpler heursitic
 	of visiblilty from a singe internal "focal" point is currently used.
-	calc_ocent() attempts to choose a point with best visibilit of
+	calc_ocent() attempts to choose a point with best visibility of
 	all the gamut surfaces, since any self-shadowing results in
 	gamut surface holes.
 
 	Improvements would be to create per-axis mappings (and separate
 	the shadow vertex locations from the real ones) to re-shape
-	the gamut into a squere as much as possible.
+	the gamut into a square as much as possible.
 	Multiple external focal points could be used,
     a vertex being shadowed only when it can't be "seen" by any
 	external focal point. It's hard to figure how to make the latter
@@ -8408,12 +8467,19 @@ static void plot_fxcell_surface(rspl *s, int dofclabels, int dobxcells, int dowa
 	accelleration structure for shadow testing (BSP tree ??),
 	and build the gamut surface incrementally from existing
 	furthest points.
+
+	(Have loop re-orderings been exausted ? i.e. can overlap
+	 triangle processing "Do a first pass for each test vertex,
+	 testing against just the triangles that are associated with
+	 it's triangle" be used for main shadowing testing ?)
 */
 
 /* Struct to slice locus points */
 struct _slpoint {
 	double v[MXRO];			/* Point location */
-	double rad;				/* Distance from ccent */
+	double wrad;			/* Weighted radius */
+	double rad;				/* Distance from ccent along slice */
+	double minrad;			/* Minimum distance from ccent */
 
 	double cvec[MXRO];		/* Vector from this point to ccent */
 	double len;				/* Length of segment, -1 if no good */
@@ -8426,12 +8492,15 @@ struct _ocenctx {
 	int ares;					/* angle resolution */ 
 	slpoint *p[MXRO];			/* Slice locus points */
 	double ccent[MXRO];			/* Construction center point */
+	double ret;					/* return value */
+	int oog;					/* flag set if center is out of gamut */
 	int debug;
 }; typedef struct _ocenctx ocenctx;
 
 /* Given a set of slice locus points and a proposed center point, */
-/* compute the weighted averag of the orthogonality of the point */
+/* compute the weighted average of the orthogonality of the point */
 /* to each locus line segment. (smaller is better) */
+/* (Used for optimizing the focal/center point.) */
 static double aorthog(void *_ctx, double *cent) {
 	ocenctx *ctx = (ocenctx *) _ctx;
 	rspl *s = ctx->s;
@@ -8441,6 +8510,8 @@ static double aorthog(void *_ctx, double *cent) {
 	double ang, aang = 0.0;
 	int naang = 0;
 	
+	ctx->oog = 0;
+
 	if (ctx->debug) printf("aorthog called with cent %s\n",debPdv(fdi,cent));
 
 	for (ff = 0; ff < fdi; ff++) {
@@ -8480,7 +8551,10 @@ static double aorthog(void *_ctx, double *cent) {
 			
 			/* Normalized difference in distance over length */
 			/* Compute dot product of cv and segment vector */
+			/* (ang range 0.0 .. 1.0 */
 			ang = fabs(trad - nrad)/ctx->p[ff][aa].len;
+			if (ang > 1.0)
+				ang = 1.0;
 
 			if (ctx->debug) printf("  aa %d: trad %f nrad %f, diff %f, len %f, ang %f\n",aa,trad,nrad,fabs(trad - nrad),ctx->p[ff][aa].len,ang);
 
@@ -8493,14 +8567,11 @@ static double aorthog(void *_ctx, double *cent) {
 
 			if (dot < 0.0) {
 				if (ctx->debug) printf("  dot is %f\n",dot);
-				ang = 50.0;			/* Big value */
+				ang = 50.0;				/* Big value */
+				ctx->oog = 1;
 			} else {
 				ang = pow(ang, 50.0);	/* Weight high angles */
 			}
-
-			if (ang > aang)
-				aang = ang;
-
 			aang += ang;
 			naang++;
 		}
@@ -8509,6 +8580,8 @@ static double aorthog(void *_ctx, double *cent) {
 	
 	if (ctx->debug) printf(" returning %f\n",aang);
 
+	ctx->ret = aang;
+
 	return aang;
 }
 
@@ -8518,6 +8591,7 @@ static double aorthog(void *_ctx, double *cent) {
 /* surface of the gamut. */
 static void calc_ocent(rspl *s) {
 	int i, j, aa, mm;
+	int e, ee, di = s->di;
 	int f, ff, fdi = s->fdi;
 	int rgres = s->rev.res;			/* number of bwd cells */
 	double minmax[2][MXRO][MXRO];	/* Range min/max points for each axis */
@@ -8527,14 +8601,15 @@ static void calc_ocent(rspl *s) {
 	double ss[MXRO];
 	ocenctx ctx;					/* Context */
 	double atanscale;
+	int ici, nici;
 
 	/* Scan the forward array for the min and max points of each axis */
 	for (f = 0; f < fdi; f++) {
-		minmax[0][f][f] = 1e30;
-		minmax[1][f][f] = -1e30;
+		minmax[0][f][f] = 1e200;
+		minmax[1][f][f] = -1e200;
 	}
 
-	/* Scan the Grid for min/max values */
+	/* Scan the fwd Grid for min/max values */
 	for (gp = s->g.a, ep = s->g.a + s->g.no * s->g.pss; gp < ep; gp += s->g.pss) {
 		for (ff = 0; ff < fdi; ff++) {
 			if (minmax[0][ff][ff] > gp[ff]) {
@@ -8550,8 +8625,8 @@ static void calc_ocent(rspl *s) {
 
 	if (fdi == 1) {
 		for (f = 0; f < fdi; f++)
-			s->rev.ocent[f] = 0.5 * (minmax[0][0][f] + minmax[1][0][f]);
-		DBG(("calc_ocent: got ocent = %s\n",debPdv(fdi,s->rev.ocent)));
+			s->rev.ocent[f] = 0.5 * (minmax[0][f][f] + minmax[1][f][f]);
+		DBG(("calc_ocent: got 1d ocent = %s\n",debPdv(fdi,s->rev.ocent)));
 		return;
 	}
 
@@ -8569,7 +8644,8 @@ static void calc_ocent(rspl *s) {
 	for (f = 0; f < fdi; f++)
 		s->rev.ocent[f] = ctx.ccent[f] /= ((fdi-1) * 2.0);
 
-	DBG(("calc_ocent: ccent = %s\n",debPdv(fdi,ctx.ccent)));
+	DBG(("calc_ocent: initial ccent = %s\n",debPdv(fdi,ctx.ccent)));
+//printf("calc_ocent: initial ccent = %s\n",debPdv(fdi,ctx.ccent));
 
 	/* If it's all to hard ... */
 	if (fdi != 3) {
@@ -8581,11 +8657,17 @@ static void calc_ocent(rspl *s) {
 		midix[f] = (int)((ctx.ccent[f] - s->rev.gl[f])/s->rev.gw[f] + 0.5);
 		mid[f] = (midix[f]+0.5) * s->rev.gw[f] + s->rev.gl[f];
 	}
+//printf("calc_ocent: mid point = %s\n",debPdv(fdi,mid));
 
 	/* Array for each slice values at angle (+ repeat at end) */
 	ctx.debug = 0;
 	ctx.s = s;
-	ctx.ares = 20;
+	ctx.ares = (rgres + 1) & ~1;	/* Make even so that there is an opposite angle */
+	if (ctx.ares < 6)
+		ctx.ares = 6;
+	else if (ctx.ares > 20)
+		ctx.ares = 20;
+//printf("  ocent ares %d\n",ctx.ares);
 	atanscale = ctx.ares/(2.0 * DBL_PI);
 	for (ff = 0; ff < fdi; ff++) {
 		if ((ctx.p[ff] = (slpoint *)rev_calloc(s, ctx.ares+1,sizeof(slpoint))) == NULL)
@@ -8593,96 +8675,201 @@ static void calc_ocent(rspl *s) {
 		INCSZ(s, (ctx.ares+1) * sizeof(slpoint));
 	}
 
-//printf("~1 locating center point\n");
+	/* Use 5 passes to locate a more reliable initial center point */
+	for (nici = 10, ici = 0; ici < nici; ici++) {
+//printf("  locating center point iter %d\n",ici);
 
-	/* Set initial radius values */
-	for (aa = 0; aa < ctx.ares; aa++) {
-		for (ff = 0; ff < fdi; ff++)
-			ctx.p[ff][aa].rad = -1.0;
-	}
-			
-	/* Take three slices through the rev[] array, plotting */
-	/* the maximum circumference for the slice */
+		/* Set initial radius values */
+		for (ff = 0; ff < fdi; ff++) {
+			for (aa = 0; aa < ctx.ares; aa++) {
+				ctx.p[ff][aa].wrad = -1.0;
+				ctx.p[ff][aa].rad = -1.0;
+				ctx.p[ff][aa].minrad = 1e38;
+			}
+		}
+				
+		/* Take three slices through the rev[] array, plotting */
+		/* the maximum circumference for the slice */
 
-	/* For the axis we're slicing */
-	for (ff = 0; ff < fdi; ff++) {
-		DCOUNT(cc, MXRO, 2, 0, 0, rgres);	/* Counter through bwd cells */
-		double vv[MXRO];
-		int aa;
-		
+		/* For the axis we're slicing */
+		for (ff = 0; ff < fdi; ff++) {
+			FCOUNT(cc, MXRO, 3);	/* Counter through bwd cells */
+			int start[3], endp1[3];
+			double vv[MXRO];
+			int aa;
+			
 //printf(" slice axis %d\n",ff);
 
-		/* Scan this slice of rev[] */
-		DC_INIT(cc);
-		while (!DC_DONE(cc)) {
-			int ix;
-			int slix[2];		/* Indexes in slice direction */
-			int *rp;
-			int co[MXRO];
-
-			/* Compute bx index */
-			ix = 0;
-			for (j = f = 0; f < fdi; f++) {
-				if (f == ff)
-					co[f] = midix[f];
-				else {
-					slix[j] = f;
-					co[f] = cc[j++];
+			/* Setup "fat" slice range */
+			for (f = 0; f < fdi; f++) {
+				if (f == ff) {
+					start[f] = midix[f]-1;
+					if (start[f] < 0)
+						start[f] = 0;
+					endp1[f] = midix[f]+2;
+					if (endp1[f] > rgres)
+						endp1[f] = rgres;
+				} else {
+					start[f] = 0;
+					endp1[f] = rgres;
 				}
-				ix += co[f] * s->rev.coi[f];
 			}
-//printf("  bx %d, %d ix %d\n",co[0],co[1],co[2],ix);
+			FRECONFA(cc, start, endp1);
+
+//printf("  slice range %d - %d, %d - %d, %d - %d\n", start[0], endp1[0]-1, start[1], endp1[1]-1, start[2], endp1[2]-1);
+
+			/* Scan this 3 thick, 2D slice of rev[] */
+			FC_INIT(cc);
+			while (!FC_DONE(cc)) {
+				int ix;
+				int slix[MXRO];		/* Indexes in slice direction + orthogonal */
+				int *rp;
+
+				/* Compute bx index */
+				ix = 0;
+				for (j = f = 0; f < fdi; f++) {
+					ix += cc[f] * s->rev.coi[f];
+					if (f != ff)
+						slix[j++] = f;
+				}
+				slix[j++] = ff;
+//printf("  bx %d, %d ix %d\n",cc[0],cc[1],cc[2],ix);
 
-			if (s->rev.rev[ix] == NULL) {
+				if (s->rev.rev[ix] == NULL) {
 //printf("  rev is empty\n");
-				goto next_bx;
-			}
+					goto next_bx;
+				}
 
-			/* For all the vertex values in bx rev[] */
-			for (rp = s->rev.rev[ix]+3; *rp != -1; rp++) {
-				float *fcb = s->g.a + *rp * s->g.pss;
-				double x, y, rad, ang;
+				/* For all the cubes bx rev[] */
+				for (rp = s->rev.rev[ix]+3; *rp != -1; rp++) {
 
-//printf("   vtx %d\n",*rp);
-				/* Ignore over ink limit values */
-				if (s->limiten && fcb[-1] > s->limitv)
-					continue;
-		
-				/* Compute radius and normalize */
-				x = fcb[slix[0]] - ctx.ccent[slix[0]];
-				y = fcb[slix[1]] - ctx.ccent[slix[1]];
-				rad = sqrt(x * x + y * y);
-				if (rad < EPS)
-					continue;
+					/* For each vertx of this cube */
+					for (ee = 0; ee < (1<<di); ee++) {
+						int vix = *rp + s->g.hi[ee];
+						float *gp = s->g.a + vix * s->g.pss;	/* Pointer to float of fwd vertex */
+						double fcb[MXRO];
+						double x, y, z, wrad, rad, ang;
 
-				/* Quantized angle this point is at */
-				ang = atanscale * atan2(y, x);
-				aa = (int)floor(ang);
-				if (aa < 0)
-					aa += ctx.ares;
-				if (aa >= ctx.ares) 
-					aa -= ctx.ares;
+						/* Don't add over ink limit vertexes */
+						if (s->limiten && gp[-1] > s->limitv) {
+							continue;
+						}
+
+						for (f = 0; f < fdi; f++)
+							fcb[f] = gp[f]; 
+
+						/* (Don't) Convert output values to log values */
+						/* logcomp(s, fcb, fcb, ctx.ccent); */
+
+						/* Compute 2D radius and normalize */
+						x = fcb[slix[0]] - ctx.ccent[slix[0]];
+						y = fcb[slix[1]] - ctx.ccent[slix[1]];
+						z = fcb[slix[2]] - ctx.ccent[slix[2]];
+						/* wrad is "elipsoid" weighted radius in slice */
+						wrad = x * x + y * y - 1.5 * z * z;
+						wrad = sqrt(wrad < 0.0 ? 0.0 : wrad);
+						rad = sqrt(x * x + y * y);
+						if (rad < EPS || wrad < EPS)
+							continue;
+
+						/* Quantized angle this point is at */
+						ang = atanscale * atan2(y, x);
+						aa = (int)floor(ang);
+						if (aa < 0)
+							aa += ctx.ares;
+						if (aa >= ctx.ares) 
+							aa -= ctx.ares;
 
 //printf("   slice %d vtx %f %f %f rad %f, ang %f aa %d\n", ff, fcb[0], fcb[1], fcb[2], rad, ang, aa);
 
-				if (rad > ctx.p[ff][aa].rad) {
-					ctx.p[ff][aa].rad = rad;
+						if (wrad > ctx.p[ff][aa].wrad) {
+							ctx.p[ff][aa].wrad = wrad;
+							ctx.p[ff][aa].rad = rad;
 
-					for (f = 0; f < fdi; f++) {
-						ctx.p[ff][aa].v[f] = fcb[f];
+							/* Copy far point */
+							for (f = 0; f < fdi; f++)
+								ctx.p[ff][aa].v[f] = fcb[f];
+
+							/* (don't) Flatten the points to lie on the notional center */
+							/* ctx.p[ff][aa].v[ff] = ctx.ccent[ff]; */
+						}
+
+						/* Track min in case ccent is not within slice */
+						if (rad < ctx.p[ff][aa].minrad) {
+							ctx.p[ff][aa].minrad = rad;
+						}
 					}
-					/* Flatten the points to lie on the notional center */
-					ctx.p[ff][aa].v[ff] = ctx.ccent[ff]; 
 				}
+			  next_bx:;
+				FC_INC(cc);
 			}
-		  next_bx:;
-			DC_INC(cc);
+
+			/* Repeat first in extra at end */
+			ctx.p[ff][ctx.ares] = ctx.p[ff][0];		/* Structure copy */
 		}
 
-		/* Repeat first in extra at end */
-		ctx.p[ff][ctx.ares] = ctx.p[ff][0];		/* Structure copy */
+		/* Check if center point is within slice by looking for empty entries. */
+		{
+			double ccvec[MXRO];		/* Center correction vector */
+			double ccount = 0.0;
+			
+			for (f = 0; f < fdi; f++)
+				ccvec[f] = 0.0;
+
+			for (ff = 0; ff < fdi; ff++) {
+				for (aa = 0; aa < ctx.ares; aa++) {
+
+//printf("   slice %d aa %d, vtx %s rad %f, irad %f\n", ff, aa, debPdv(fdi,ctx.p[ff][aa].v), ctx.p[ff][aa].rad, ctx.p[ff][aa].minrad);
+
+					/* Either the grid is very sparse, or our center */
+					/* is outside */
+					if (ctx.p[ff][aa].rad < 0.0) {
+						int oaa = aa + (ctx.ares/2);
+						if (oaa >= ctx.ares) 
+							oaa -= ctx.ares;
+
+//printf("   oaa %d, vtx %s rad %f, irad %f\n", oaa, debPdv(fdi,ctx.p[ff][oaa].v), ctx.p[ff][oaa].rad, ctx.p[ff][oaa].minrad);
+
+						/* If oposite side has an entry */
+						if (ctx.p[ff][oaa].rad > 0.0) {
+							double cor[MXRO];
+							double prop = (3.0 * ctx.p[ff][oaa].minrad
+							                   + ctx.p[ff][oaa].rad)/4.0;
+
+							prop /= ctx.p[ff][oaa].rad;		/* Proportion of distance to v[] */
+
+							for (f = 0; f < fdi; f++)
+								cor[f] = prop * (ctx.p[ff][oaa].v[f] - ctx.ccent[f]);
+//printf("   prop %f, corr %s\n",prop,debPdv(fdi,cor));
+							for (f = 0; f < fdi; f++)
+								ccvec[f] += cor[f];
+							ccount++;
+						}
+					}
+				}
+			}
+
+//printf("ccount %f\n",ccount);
+			if (ccount > 0.0) {	/* Make adjustment */
+				if (ici < (nici-1)) {
+					for (f = 0; f < fdi; f++)
+						ccvec[f] /= ccount;
+//printf("Corecting center by %s\n",debPdv(fdi,ccvec));
+					for (f = 0; f < fdi; f++)
+						ctx.ccent[f] += ccvec[f];
+//printf("cceny now %s\n",debPdv(fdi,ctx.ccent));
+				} else {		/* Last round and correction needed */
+					if (0.0 && s->verbose)
+						fprintf(stdout, "%cFailed to locate aprox. gamut center\n",cr_char);
+				}
+			} else {
+				break;		/* We're done */
+			}
+		}
 	}
 
+//printf("calc_ocent: refined ccent = %s\n",debPdv(fdi,ctx.ccent));
+
 	/* Pre-compute point to point info to speed optimization */
 	for (ff = 0; ff < fdi; ff++) {
 		for (aa = 0; aa < ctx.ares; aa++) {
@@ -8706,7 +8893,9 @@ static void calc_ocent(rspl *s) {
 	/* slice segment. This should maximize visibility of the inner of the */
 	/* gamut surface, for shadow testing. */
 	for (f = 0; f < fdi; f++)
-		ss[f] = 5.0;
+		ss[f] = fabs(0.1 * (minmax[1][f][f] - minmax[0][f][f]));
+
+//ctx.debug = 1;
 
 	/* return 0 on sucess, 1 on failure due to excessive itterations */
 	if (powell(NULL, fdi, s->rev.ocent, ss, 1e-3, 500, aorthog, (void *)&ctx, NULL, NULL)) {
@@ -8715,8 +8904,16 @@ static void calc_ocent(rspl *s) {
 			s->rev.ocent[f] = ctx.ccent[f];
 	}
 
-//	ctx.debug = 1;
-//	printf("Final angle = %f\n", aorthog(&ctx, ctx.ccent));
+//ctx.debug = 1;
+
+	/* Check result */
+	aorthog(&ctx, ctx.ccent);
+
+	/* Hmm. This isn't very reliable in detecting failure. */
+	if (ctx.oog)
+		printf("calc_ocent failed to return in-gamut focal point!\n");
+			
+//printf("Final angle = %f\n", ctx.ret);
 
 #ifdef REVVRML		/* Plotting routine declarations */
 	/* Diagnostic - dump the gamut slice locii */
@@ -8729,7 +8926,7 @@ static void calc_ocent(rspl *s) {
 		double blue[3] = { 0.1, 0.1, 0.8 };
 		double *rgb[3] = { red, green, blue };
 
-		wrl = new_vrml("section", 0, vrml_lab);
+		wrl = new_vrml("section", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
 		wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
 
 		/* Show vertex labels */
@@ -8772,7 +8969,384 @@ static void calc_ocent(rspl *s) {
 		DECSZ(s, (ctx.ares+1) * sizeof(slpoint));
 	}
 
-	DBG(("calc_ocent: got ocent = %s\n",debPdv(fdi,s->rev.ocent)));
+	DBG(("calc_ocent: final ocent = %s\n",debPdv(fdi,s->rev.ocent)));
+//printf("calc_ocent: final ocent = %s\n",debPdv(fdi,s->rev.ocent));
+}
+
+/* Create gamut surface linearization (surflin) transform. */
+/* This is used by logcomp() to try and straighten out the */
+/* device response so that the ocent is "visible" from */
+/* any point on the surface. */ 
+/* (We assume we are called at the correct point when bx->status == bx_uninit) */
+static int calc_surflin(
+rspl *s,
+vtxcache *vc,		/* Vertexes */
+assdire *edgdir		/* Edge lookup for vertex */
+) {
+	int i, j, g;
+	int e, ee, di = s->di;
+	int f, ff, fdi = s->fdi;
+	vtxrec *vx, *nvx;
+	int nitter, itter;
+
+	double vxv[POW2MXRI][MXRO];		/* Overal fwd interp vertex values */
+	int nvtx;
+	cow *mpoints;
+	int gres[MXDO];
+	double min[MXDO], max[MXDO];
+	double vmin[MXDO], vmax[MXDO];
+
+	for (f = 0; f < fdi; f++)
+		gres[f] = s->g.bres;		// ??
+
+//printf("calc_surflin: rspl res %d\n",gres[0]);
+	DBG(("calc_surflin: rspl res %d\n",gres[0]));
+
+//printf("~1 gres = %d %d %d\n", s->g.res[0], s->g.res[1], s->g.res[2]);
+//printf("~1 ci = %d %d %d\n", s->g.ci[0], s->g.ci[1], s->g.ci[2]);
+
+	/* Lookup interpolation cube corners */
+	for (ee = 0; ee < (1<<di); ee++) {
+		int ix;
+		float *fcb;
+		for (ix = e = 0; e < di; e++) {
+			if (ee & (1<<e))
+				ix += s->g.ci[e] * (s->g.res[e] - 1);
+		}
+		fcb = s->g.a + ix * s->g.pss;
+		for (f = 0; f < fdi; f++)
+			vxv[ee][f] = fcb[f];
+//printf("~1 cube corners %d = %f %f %f\n",ee, vxv[ee][0], vxv[ee][1], vxv[ee][2]);
+	}
+
+//printf("calc_surflin: counting number of vertexes:\n");
+	/* Count the number of vertexes we may need */
+	nvtx = 0;
+	for (i = 0; i < vc->hash_size; i++) {
+		for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+			nvtx++;
+		}
+	}
+	DBG(("calc_surflin: %d mapping points\n",nvtx));
+//printf("calc_surflin: %d mapping points\n",nvtx);
+
+//printf("calc_surflin: computing goal values:\n");
+	/* Currently the vertex vl = vv = output value of vertex. */
+	/* Temporarily replace vv with the idealize (linear interp) output "goal" values. */
+	for (i = 0; i < vc->hash_size; i++) {
+		for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+			int tix;				/* Temp fwd cell index */
+			double we[MXRI];		/* Vertex input position == 1.0 - Weight */
+			double gw[POW2MXRI];	/* weight for each grid cube corner */
+			double w;
+
+			/* Compute this vertexes relative input position */
+			for (tix = vx->ix, e = 0; e < di; e++) {
+				int dix;
+				dix = tix % s->g.res[e];
+				tix /= s->g.res[e];
+				we[e] = (double)dix/(s->g.res[e]-1.0);
+			}
+			
+			/* Compute corner weights needed for interpolation */
+			gw[0] = 1.0;
+			for (e = 0, g = 1; e < di; g *= 2, e++) {
+				for (j = 0; j < g; j++) {
+					gw[g+j] = gw[j] * we[e];
+					gw[j] *= (1.0 - we[e]);
+				}
+			}
+		
+			/* Linear interpolated output values */
+			w = gw[0];
+			for (f = 0; f < fdi; f++)			/* Base of cube */
+				vx->vv[f] = w * vxv[0][f];
+		
+			for (g = 1; g < (1<<di); g++) {	/* For all other corners of cube */
+				w = gw[g];
+				for (f = 0; f < fdi; f++)
+					vx->vv[f] += w * vxv[g][f];
+			}
+		}
+	}
+
+	/* Now itteratively adjust the vl values to better match the scaled */
+	/* relative positions of the goal values. */
+
+	/* Go through them again to get every line they are part of */
+	/* (We're assuming we need exponentially more itters with finer point */
+	/* spacing ?) */
+	nitter = (int)(0.06 * s->g.bres * s->g.bres + 0.5);
+	if (nitter < 1)
+		nitter = 1;
+	for (itter = 0; itter < nitter; itter++) {
+
+		DBG(("calc_surflin: maping itter %d\n",itter));
+//printf("calc_surflin: maping itter %d\n",itter);
+
+		for (i = 0; i < vc->hash_size; i++) {
+			for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+				assdire *edg;			/* Edge table */
+				float *fp;
+				int fl;
+				double agrad, aorad;
+				int nn;
+				double scale;
+
+				fp = s->g.a + vx->ix * s->g.pss;	/* This vertex in fwd grid */
+				fl = FLV(fp);		/* Edge flags for this vertex */
+				edg = edgdir + fl;
+
+				/* For vertexes at the end of all possible edges common with this vertex, */
+				/* compute average radius */
+				agrad = aorad = 0.0;
+				for (j = 0; j < edgdir[fl].no; j++) {
+					int eix;
+
+					/* Index number of vertex other than the one we got it from */
+					if (edg->ti[j].goffs[0] != 0)
+						eix = vx->ix + edg->ti[j].goffs[0];
+					else
+						eix = vx->ix + edg->ti[j].goffs[1];
+
+					if ((nvx = get_vtxrec(vc, eix)) != NULL) {
+						double glen, olen;
+						glen = olen = 0.0;
+						for (f = 0; f < fdi; f++) {
+							double tt;
+
+							tt = vx->vv[f] - nvx->vv[f];
+							glen += tt * tt;
+
+							tt = vx->vl[f] - nvx->vl[f];
+							olen += tt * tt;
+						}
+						glen = sqrt(glen);
+						olen = sqrt(olen);
+						agrad += glen;
+						aorad += olen;
+						nn++;
+					}
+				}
+
+				if (nn == 0) {					/* Hmm. No neighors ? */
+					vx->status = vtx_del;		/* Mark it as isolated */
+					continue;
+				}
+
+				scale = aorad/agrad;	/* Local scale factor from goal to output */
+
+				/* Reset the current vertex output value based on the relative */
+				/* position of it in goal space */
+				for (f = 0; f < fdi; f++)
+					vx->vl[f] = 0.0;
+
+				nn = 0;
+				for (j = 0; j < edgdir[fl].no; j++) {
+					int eix;
+
+					/* Index number of vertex other than the one we got it from */
+					if (edg->ti[j].goffs[0] != 0)
+						eix = vx->ix + edg->ti[j].goffs[0];
+					else
+						eix = vx->ix + edg->ti[j].goffs[1];
+
+					if ((nvx = get_vtxrec(vc, eix)) != NULL) {
+						for (f = 0; f < fdi; f++)
+							vx->vl[f] += nvx->vl[f] + scale * (vx->vv[f] - nvx->vv[f]);
+						nn++;
+					}
+				}
+				for (f = 0; f < fdi; f++)
+					vx->vl[f] /= (double)nn;
+			}
+		}
+	}
+
+	DBG(("calc_surflin: creating rspl\n"));
+//printf("calc_surflin: creating rspl\n");
+
+	/* Now construct rspl setup mapping points from vertex normal output values */
+	/* to adjusted vl values */
+
+	if ((s->rev.surflin = new_rspl(RSPL_NOFLAGS, fdi, fdi)) == NULL)
+		error("calc_surflin: new_rspl failed");
+
+	nvtx++;		/* One for center point */
+
+	/* Allocate rspl setup points */
+	if ((mpoints = malloc(sizeof(cow) * 2 * nvtx)) == NULL)
+//	if ((mpoints = malloc(sizeof(cow) * nvtx)) == NULL)
+		error("calc_surflin: malloc of %d rspl setup points failed",nvtx);
+	
+	nvtx = 0;
+
+#ifndef NEVER
+	/* Center point */
+	for (f = 0; f < fdi; f++) {
+		mpoints[nvtx].p[f] = s->rev.ocent[f]; 
+		mpoints[nvtx].v[f] = s->rev.ocent[f]; 
+	}
+	mpoints[nvtx].w = 10.0;
+	nvtx++;
+#endif
+
+	/* Set the surface mapping points and restore vertexes values */
+	for (i = 0; i < vc->hash_size; i++) {
+		for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+			float *fcb;
+
+			fcb = s->g.a + vx->ix * s->g.pss;	/* This vertex in fwd grid */
+
+			/* Actual output value as source of mapping */
+			for (f = 0; f < fdi; f++)
+				mpoints[nvtx].p[f] = vx->vv[f] = fcb[f];
+
+			if (vx->status != vtx_norm) {	/* Skip isolated values */
+				vx->status = vtx_norm;		/* Restore vtx contents */
+				for (f = 0; f < fdi; f++)
+					vx->vl[f] = vx->vv[f];
+				continue;
+			}
+
+			for (f = 0; f < fdi; f++)
+				mpoints[nvtx].v[f] = vx->vl[f];
+			mpoints[nvtx].w = 1.0;
+			nvtx++;
+
+			vx->status = vtx_norm;		/* Restore vtx contents */
+			for (f = 0; f < fdi; f++)
+				vx->vl[f] = vx->vv[f];
+			
+//printf("Map[%d] %f %f %f -> %f %f %f\n", nvtx-1, mpoints[nvtx-1].p[0], mpoints[nvtx-1].p[1], mpoints[nvtx-1].p[2], mpoints[nvtx-1].v[0], mpoints[nvtx-1].v[1], mpoints[nvtx-1].v[2]);
+
+#ifndef NEVER
+			/* Add intermediate "fixed" point */
+			for (f = 0; f < fdi; f++) {
+				mpoints[nvtx].p[f] = mpoints[nvtx].v[f]
+				                  = 0.5 * mpoints[nvtx-1].p[f] + 0.5 * s->rev.ocent[f];
+			}
+			mpoints[nvtx].w = 0.5;
+			nvtx++;
+#endif
+		}
+	}
+
+	for (f = 0; f < fdi; f++) {
+		min[f] = 1e38;
+		max[f] = -1e38;
+		vmin[f] = 1e38;
+		vmax[f] = -1e38;
+	}
+	for (i = 0; i < nvtx; i++) {
+
+#ifdef NEVER
+		/* Blend with original values */
+		for (f = 0; f < fdi; f++)
+			mpoints[i].v[f] = 0.5 * mpoints[i].p[f] + 0.5 * mpoints[i].v[f];
+#endif
+
+		for (f = 0; f < fdi; f++) {
+			if (mpoints[i].p[f] < min[f])
+				min[f] = mpoints[i].p[f];
+			if (mpoints[i].p[f] > max[f])
+				max[f] = mpoints[i].p[f];
+
+			if (mpoints[i].v[f] < vmin[f])
+				vmin[f] = mpoints[i].v[f];
+			if (mpoints[i].v[f] > vmax[f])
+				vmax[f] = mpoints[i].v[f];
+		}
+	}
+
+#ifdef REVVRML	/* Plot mapping vectors red->green */
+	{
+		vrml *wrl;
+		double red[3] = { 1.0, 0.0, 0.0 };
+		double green[3] = { 0.0, 1.0, 0.0 };
+
+		wrl = new_vrml("suflinvecss", 0, vrml_lab);
+		wrl->start_line_set(wrl, 0);
+
+		for (i = 0; i < nvtx; i++) {
+			wrl->add_col_vertex(wrl, 0, mpoints[i].p, red);
+			wrl->add_col_vertex(wrl, 0, mpoints[i].v, green);
+			
+		}
+
+		wrl->make_lines(wrl, 0, 2);
+		wrl->del(wrl);
+	}
+#endif
+
+	DBG(("calc_surflin: mapping points set, about to creat rspl:\n"));
+//printf("calc_surflin: mapping points set, about to creat rspl:\n");
+
+	/* Fit the rspl */
+	s->rev.surflin->fit_rspl_w(s->rev.surflin, RSPL_NOFLAGS, mpoints, nvtx,  
+	                           min, max, gres, vmin, vmax, 4.0, NULL, NULL);   
+
+	DBG(("calc_surflin: mapping created\n"));
+//printf("calc_surflin: mapping created\n");
+
+#ifdef NEVER
+	{
+		double de;
+		co p;
+		extern double icmNorm33(double *, double *);
+
+		/* Check fit */
+		de = 0.0;
+
+		for (i = 0; i < nvtx; i++) {
+			for (f = 0; f < fdi; f++)
+				p.p[f] = mpoints[i].p[f];
+
+			s->rev.surflin->interp(s->rev.surflin, &p);
+	
+			de += icmNorm33(mpoints[i].v, p.v);
+		}
+		de = de/(double)nvtx;
+
+		printf("Avg fit error = %f\n",de);
+	}
+#endif
+	
+	free(mpoints);
+
+	/* Lookup ocent mapping offset */
+	{
+		co p;
+
+		for (f = 0; f < fdi; f++)
+			p.p[f] = s->rev.ocent[f];
+		s->rev.surflin->interp(s->rev.surflin, &p);
+//printf("opoint mapping %f %f %f -> %f %f %f\n", p.p[0], p.p[1], p.p[2], p.v[0], p.v[1], p.v[2]);
+		for (f = 0; f < fdi; f++)
+			s->rev.linoff[f] = p.v[f] - s->rev.ocent[f];
+	}
+	
+
+#ifndef NEVER
+	/* Put the transform into use */
+	s->rev.surflin_en = 1;
+
+	/* Transform all the vertexes */
+	for (i = 0; i < vc->hash_size; i++) {
+		for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+			logcomp(s, vx->vl, vx->vv, s->rev.ocent);
+
+			/* Compute distance to overall center point squared */
+			vx->dist = 0.0;
+			for (f = 0; f < fdi; f++) {
+				double tt = vx->vl[f] - s->rev.ocent[f]; 
+				vx->dist += tt * tt;
+			}
+		}
+	}
+#endif
+
+	return 0;
 }
 
 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
@@ -8800,6 +9374,7 @@ rspl *s
 								/* Note that bit 1 can be set for cells that are not */
 								/* to be explored because they are in the gamut interior, */
 								/* and because they have already been added to the seedlist. */
+	int pass;					/* Construction pass */
 	float *gp;					/* Pointer to fwd grid points */
 
 	DCOUNT(gg, MXRO, fdi, 0, 0, rgres);	/* Track the prime seed coordinate */
@@ -8835,9 +9410,24 @@ rspl *s
 	int nnmxrevcelldepth = 0;			/* Maximum nnrev[] list lengths */
 	int nnrevshare = 0;					/* Sum of nnrev[] list reference counts */
 #endif
+	datao rgmin, rgmax;
 
 	DBG(("init_revaccell called, di = %d, fdi = %d, mgres = %d\n",di,fdi,(int)s->g.mres));
 
+	/* To help VRML diagnostics, make a guess as to whether the output */
+	/* space is XYZ like, or L*a*b* like */
+
+	s->get_out_range(s, rgmin, rgmax);	/* overall output min/max */
+
+	if (fdi >= 3
+	 && rgmin[0] >= -1.0 && rgmax[0] < 3.0
+	 && rgmin[1] >= -1.0 && rgmax[1] < 3.0
+	 && rgmin[2] >= -1.0 && rgmax[2] < 3.0) {
+		s->rev.probxyz = 1;
+		if (s->verbose)
+			fprintf(stdout, "%cLooks like an XYZ space\n",cr_char);
+	}
+
 	if (fdi > 1 && s->verbose)
 		fprintf(stdout, "%cInitializing nnrev arrays...\n",cr_char);
 
@@ -8885,7 +9475,7 @@ rspl *s
 	/*
 	 * The rev[] and nnrev[] grids contain pointers to lists of grid cube base indexes.
 	 * If the pointer is NULL, then there are no base indexes in that list.
-	 * A non NULL list uses element [0] to indicate the alocation size of the list,
+	 * A non NULL list uses element [0] to indicate the allocation size of the list,
 	 * [1] contains the index of the next free location, [2] contains the reference
      * count (lists may be shared), the list starts at [3]. The last entry is marked with -1.
 	 */
@@ -8897,6 +9487,8 @@ rspl *s
 		double iv[MXDI];				/* Input value corresponding to grid */
 
 		DBG(("Looking up fwd vertex ink limit values\n"));
+//printf("Looking up fwd vertex ink limit values\n");
+//printf("s->limitv = %f\n",s->limitv);
 		/* Calling the limit function for each fwd vertex could be bad */
 		/* if the limit function is slow. Maybe an octree type algorithm */
 		/* could be used if this is a problem ? */
@@ -8906,7 +9498,9 @@ rspl *s
 				for (e = 0; e < di; e++)
 					iv[e] = s->g.l[e] + gc[e] * s->g.w[e];  /* Input sample values */
 				gp[-1] = (float)(INKSCALE * s->limitf(s->lcntx, iv));
+//printf("~1 set ix %d limitv to %f\n",i,gp[-1]);
 			}
+//else printf("~1 ix %d limitv is %f\n",i,gp[-1]);
 			EC_INC(gc);
 		}
 		s->g.limitv_cached = 1;
@@ -8958,6 +9552,7 @@ rspl *s
 		uil = oil = 0;
 		for (f = 0; f < fdi; f++)	/* Init output min/max */
 			min[f] = max[f] = gp[f];
+
 		if (!s->limiten || gp[-1] <= s->limitv)
 			uil = 1;
 		else
@@ -9065,7 +9660,7 @@ rspl *s
 		}	/* Next reverse grid point in intersecting cube */
 	}	/* Next base grid point */
 
-	DBG(("We skipped %d cells that were over the limit\n",nskcells));
+	DBG(("We skipped %d/%d cells that were over the limit\n",nskcells,gno));
 
 #ifdef CHECK_NNLU
 	if (fdi > 1) {
@@ -9180,7 +9775,6 @@ rspl *s
 				cc[f]  += (ff & 1) ? 1 : -1;
 				vflagp += (ff & 1) ? s->rev.coi[f] : -s->rev.coi[f];
 
-
 				/* Out of bounds or empty */
 				if (cc[f] < 0 || cc[f] >= rgres || ((*vflagp & 0xf) == 0)) {
 					vflag[i] = (vflag[i] & ~0xf) | 2;	/* Convert this one to empty surface cell */
@@ -9256,7 +9850,7 @@ rspl *s
 			}
 		}
 
-	} else {
+	} else {		/* di >= 3 gamut surface finding */
 
 #ifdef REVVRML
 		/* Plot the initial surface bxcells & their fwd cells. */
@@ -9289,7 +9883,7 @@ rspl *s
 		/* - - - - - - - - - - - - - - */
 		/* fill, thin and add, until */
 		/* there is no more work to do. */
-		for (;;) {
+		for (pass = 0;; pass++) {
 			int phase;
 			int morevtxadded = 0;
 #if defined(REVTABLESTATS) || defined(DEBUG)
@@ -9301,7 +9895,7 @@ rspl *s
 
 			/* For each surface bxcell, convert the corresponding */
 			/* rev[] fwd cubes into vertices. */
-			/* (Musk keep bxcells even if none of their verticies */
+			/* (Must keep bxcells even if none of their verticies */
 			/*  are physically in them, so that those verticies get thinned. */
 			/*  could only remove them if vertex was not in any surface cell ?) */
 			for (pbx = &s->rev.surflist, bx = *pbx; bx != NULL; bx = nbx) {
@@ -9368,10 +9962,10 @@ rspl *s
 				/* Expand a bxcell's shadow testing group values based on it's vertex list */
 				/* so that shadow testing works correctly for vertexes that don't */
 				/* actually lie within the bxcell. (Note that in fact the triangle */
-				/* testing creates triangles that are mode of vertexes that may not */
-				/* be in this bx's list, so the shadow size doesn't accturatly  reprsent */
+				/* testing creates triangles that are made of vertexes that may not */
+				/* be in this bx's list, so the shadow size doesn't accuratly  reprsent */
 				/* the possible shadow area. It's not clear what consequences this has, */
-				/* if any. If we extanded the group to cover this, we would need to have ) */
+				/* if any. If we extanded the group to cover this, we would need to have */
 				/* two groups, a shadower group including those vertexes, and a shadowee */
 				/* goup for just those vertexes that are part of the bx. */
 				extend_bxcell_shadow_group(s, &vc, bx);
@@ -9380,6 +9974,19 @@ rspl *s
 				morevtxadded = 1;
 			}
 	
+			/* Compute transform rspl that helps "unfold" any regions of the surface */
+			/* that overlap from the perspective of ocent, to try and avoid gaps in */
+			/* the final gamut surface. Existing vtxrec are converted to have vl */
+			/* in the unfolded space. */
+			if (pass == 0 /* && function flag set */) {
+#ifndef EN_UNTWIST		/* Control using an environment variable */
+				if (getenv("ARGYLL_UNTWIST_GAMUT_SURFACE") != NULL)
+#endif
+				{
+					calc_surflin(s, &vc, edgdir);
+				}
+			}
+
 			DBG(("thinning surface vertex lists and converting to cells\n"));
 
 			/* (Sorting bxcells doesn't seem to make any performace difference.) */
@@ -9440,8 +10047,8 @@ rspl *s
 								/* If any of bx is further from nbx and their bounding */
 								/* cylinders overlap in perspective from rev.ocenter, */
 								/* assume nbx is a shadow */
-								if (shadow_group_group(s, s->rev.ocent, bx->g.bcent, bx->cc, bx->dw,    
-									                                    nbx->g.bcent, nbx->cc, nbx->dw)) {
+								if (shadow_group_group(s, s->rev.ocent, bx->g.bcent, bx->cc,
+									                 bx->dw, nbx->g.bcent, nbx->cc, nbx->dw)) {
 									nbx->wlist = bx->wlist;
 									bx->wlist = nbx;
 //printf("~1 adding shadow nnrev[%d] from xlist\n",nbx->ix);
@@ -9511,8 +10118,8 @@ rspl *s
 							/* marked deleted ??) */
 							if (
 //						    vx->status == vtx_norm &&
-							    shadow_group_vertex(s, s->rev.ocent, bx->g.bcent, bx->cc, bx->dw,
-								                                                          vx->v)) { 
+							    shadow_group_vertex(s, s->rev.ocent, bx->g.bcent, bx->cc,
+								                                          bx->dw, vx->vl)) { 
 									add_vtxrec_list(&vc, vx, 0);	/* Add if not deleted */
 //printf(" Added ix %d from bx %d\n",vx->ix,nbx->ix);
 							}
@@ -9529,7 +10136,8 @@ rspl *s
 							error("Failed to find vertex %s in cache",*rp);
 
 						if (vx->status == vtx_norm &&
-						    shadow_group_vertex(s, s->rev.ocent, bx->g.bcent, bx->cc, bx->dw, vx->v)) { 
+						    shadow_group_vertex(s, s->rev.ocent, bx->g.bcent, bx->cc,
+							                                          bx->dw, vx->vl)) { 
 								add_vtxrec_list(&vc, vx, 1);	/* Add if not hidden/deleted */
 						}
 					}
@@ -9543,7 +10151,7 @@ rspl *s
 					/* For vertexes of this bxcell and shadowers, */
 					/* in order from largst to smallest distance from center. */
 					for (vx = vc.vtxlist; vx != NULL; vx = vx->tlist) {
-						float *vp;				/* Vertex being tested */
+						float *fcb;				/* Vertex being tested */
 						int fl;
 						assdire *tri;			/* Triangle table */
 
@@ -9556,8 +10164,8 @@ rspl *s
 
 //printf("~1 doing vertex %d at %s dist %f\n",vx->ix, debPdv(fdi,vx->v), sqrt(vx->dist));
 
-						vp = s->g.a + vx->ix * s->g.pss;	/* This vertex in fwd grid */
-						fl = FLV(vp);		/* Edge flags for this vertex */
+						fcb = s->g.a + vx->ix * s->g.pss;	/* This vertex in fwd grid */
+						fl = FLV(fcb);		/* Edge flags for this vertex */
 
 						tri = tridir + fl;
 //printf("~1 fl %d = 0o%o, no triangles %d\n",fl, fl, tri->no);
@@ -9597,7 +10205,7 @@ rspl *s
 
 							for (e = 0; e <= sdi; e++) {
 								for (f = 0; f < fdi; f++)
-									v[e][f] = trivx[e]->v[f];
+									v[e][f] = trivx[e]->vl[f];
 							}
 
 							/* Compute shadow group params of triangle for quick vertex test */
@@ -9633,7 +10241,7 @@ rspl *s
 
 								/* Do quick check against triangle */
 								if (!shadow_group_vertex(s,
-									     s->rev.ocent, gc, cc, dw, nvx->v)) { 
+									     s->rev.ocent, gc, cc, dw, nvx->vl)) { 
 //printf("~1 shadow group check shows no intersection\n");
 									continue;
 								}
@@ -9644,11 +10252,14 @@ rspl *s
 								shdwd = wsrv = 0;
 
 								/* Compute line delta */
-								vp = s->g.a + nvx->ix * s->g.pss;
-								for (f = 0; f < fdi; f++) {
-									pv[f] = vp[f];
+								fcb = s->g.a + nvx->ix * s->g.pss;
+								for (f = 0; f < fdi; f++)
+									pv[f] = fcb[f];
+
+								logcomp(s, pv, pv, s->rev.ocent);
+
+								for (f = 0; f < fdi; f++)
 									de[f] = pv[f] - s->rev.ocent[f];
-								}
 
 								/* Setup line cla and clb */
 								init_line_eq_imp(s, NULL, &cla, clb, s->rev.ocent, de, 0);
@@ -9773,13 +10384,13 @@ rspl *s
 						/* Check all of its neighbor vertexes, to see if */
 						/* it's safe to actually delete them. */
 						if (vx->status >= vtx_sha) {	/* vertex to delete ? */
-							float *vp;
+							float *fcb;
 							int fl;
 							assdire *edg;			/* Edge table */
 
 //printf("Checking vx %d neighbors\n",vx->ix);
-							vp = s->g.a + vx->ix * s->g.pss;	/* This vertex in fwd grid */
-							fl = FLV(vp);		/* Edge flags for this vertex */
+							fcb = s->g.a + vx->ix * s->g.pss;	/* This vertex in fwd grid */
+							fl = FLV(fcb);		/* Edge flags for this vertex */
 							edg = edgdir + fl;
 
 							/* For all possible edges that use this vertex */
@@ -9972,7 +10583,7 @@ rspl *s
 
 			/* Struct to hold test vertex locations */
 			struct _tvxrec {
-				double v[MXRO];
+				double v[MXRO];			/* Log output vertex value */
 				double dist;			/* Distance from center point squared */
 				int ix[MXRO+1];			/* Indexes of the triangle verticies */
 				int shad;				/* Test result */
@@ -10082,10 +10693,10 @@ rspl *s
 						for (j = 0; j <= sdi;  j++) {
 							if (trivx[j]->status == vtx_norm) {
 								for (f = 0; f < fdi; f++)
-									tvx->v[f] += 0.95/nntvsh * trivx[j]->v[f];
+									tvx->v[f] += 0.95/nntvsh * trivx[j]->vl[f];
 							} else {
 								for (f = 0; f < fdi; f++)
-									tvx->v[f] += 0.05/ntvsh * trivx[j]->v[f];
+									tvx->v[f] += 0.05/ntvsh * trivx[j]->vl[f];
 								trivx[j]->cross = 1;		/* For diagnostics */
 							}
 						}
@@ -10202,7 +10813,7 @@ rspl *s
 
 							for (j = 0; j <= sdi; j++) {
 								for (f = 0; f < fdi; f++)
-									v[j][f] = trivx[j]->v[f];
+									v[j][f] = trivx[j]->vl[f];
 							}
 
 							/* Compute shadow group params of triangle for quick vertex test */
@@ -10380,7 +10991,7 @@ rspl *s
 
 						for (j = 0; j <= sdi; j++) {
 							for (f = 0; f < fdi; f++)
-								v[j][f] = trivx[j]->v[f];
+								v[j][f] = trivx[j]->vl[f];
 						}
 
 						/* Compute shadow group params of triangle for quick vertex test */
@@ -10570,7 +11181,7 @@ rspl *s
 				vtxrec *vx;
 
 				if ((vx = get_vtxrec(&vc, *rp)) == NULL)
-					continue;		/* Hmm. */
+					continue;		/* Hmm. Delete it.*/
 
 				/* Keep all the un-shadowed or preserved vertexes */
 				if (vx->status == vtx_norm	
@@ -10618,9 +11229,9 @@ rspl *s
 		/* Add extra over ink limit vertexes. */
 		for (bx = s->rev.surflist; bx != NULL; bx = bx->slist) {
 			int sdi = 1;					/* sub-simplexes are edges */
-			int *crp, *rp, *nrp;
-			int ttouch;
+			int *rp;
 			vtxrec *vx, *nvx;
+			int *exlist = NULL;
 	
 			/* Add over ink limit vertexes, so that fwd cells will straddle */
 			/* the ink limit boundary. */
@@ -10630,8 +11241,8 @@ rspl *s
 			/* over ink limit verticies.) */ 
 			if (s->limiten && vflag[bx->ix] & 0x10) {
 				int *rp;
-//printf("~1 ink limitin is enabled bx %d\n", bx->ix);
 
+//printf("~1 ink limitin is enabled bx %d\n", bx->ix);
 				for (rp = bx->sl+3; *rp != -1; rp++) {
 					float *vp, *evp;
 					int fl;
@@ -10671,7 +11282,7 @@ rspl *s
 
 //printf(" Checking edge %d (%f) -> %d (%f)\n", vx->ix, vp[-1], eix, evp[-1]);
 
-						/* If over limit, add it to the list */
+						/* If over limit, add it to the expansion list */
 						if (evp[-1] > s->limitv) {
 //printf("~1 added over ink limit vertex %d\n",eix);
 
@@ -10679,10 +11290,18 @@ rspl *s
 								continue;		/* Added by another bx */
 							nvx = new_vtxrec(s, &vc, eix);
 							nvx->status = vtx_oil; 
-							add2indexlist(s, &bx->sl, eix, 0);
+							add2indexlist(s, &exlist, eix, 0);
 						}
 					}
 				}
+		
+				/* If we found over ink limit verticies, add them to our list */
+				if (exlist != NULL) {
+					for (rp = exlist+3; *rp != -1; rp++) {
+						add2indexlist(s, &bx->sl, *rp, 0);
+					}
+					free_indexlist(s, &exlist);
+				}
 			}
 		}
 
@@ -10808,6 +11427,11 @@ rspl *s
 			bx->status = bx_conv;
 		}
 
+		if (s->rev.surflin != NULL) {	/* Don't need surflin anymore */
+			s->rev.surflin->del(s->rev.surflin);
+			s->rev.surflin = NULL;
+			s->rev.surflin_en = 0;
+		}
 		if (cla != NULL)
 			free_dmatrix(cla, 0, fdi-1, 0, fdi);
 		free_trirec(s, &stc);
@@ -11911,7 +12535,7 @@ int dofwlabels			/* Plot fwd cell base indexs */
 	double grey[3] = { 0.5, 0.5, 0.5 };
 	double white[3] = { 1.0, 1.0, 1.0 };
 
-	wrl = new_vrml("raw_bxfwcells", 0, vrml_lab);
+	wrl = new_vrml("raw_bxfwcells", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
 	wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
 
 	if (dofwlabels) {
@@ -12080,7 +12704,7 @@ double xv[MXRO]			/* Intersection point */
 	double blue[3] = { 0.1, 0.1, 0.8 };
 	double yellow[3] = { 0.8, 0.8, 0.1 };
 
-	wrl = new_vrml("tri_check", 0, vrml_lab);
+	wrl = new_vrml("tri_check", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
 
 	/* Gamut center point marker */
 	wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
@@ -12166,6 +12790,8 @@ double xv[MXRO]			/* Intersection point */
 /* Main summary plot at each thinning round and at end. */
 /* Show vertex surface & optional added or deleted vertexes, */
 /* + optional bxcells. */
+#define VV vv		/* Actual surface values */
+//#define VV vl		/* Logf mapped surface values */
 static void plot_vtx_surface(
 rspl *s,
 int dovtxlabels,		/* Show vertex index numbers */
@@ -12193,9 +12819,9 @@ assdire *edgdir		/* Edge lookup for vertex */
 	bxcell *vbx;
 
 	if (dopres)
-		wrl = new_vrml("final_surface", 0, vrml_lab);
+		wrl = new_vrml("last_surface", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
 	else
-		wrl = new_vrml("thinned_surface", 0, vrml_lab);
+		wrl = new_vrml("thinned_surface", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
 	wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
 
 	if (dovtxlabels) {
@@ -12209,7 +12835,7 @@ assdire *edgdir		/* Edge lookup for vertex */
 				 || (dopres && vx->pres)
 				 || (dooil && vx->status == vtx_oil)) {
 					sprintf(index, "%d",vx->ix);
-					wrl->add_text(wrl, index, vx->v, cyan, 0.3);
+					wrl->add_text(wrl, index, vx->VV, cyan, 0.3);
 				}
 			}
 		}
@@ -12218,20 +12844,25 @@ assdire *edgdir		/* Edge lookup for vertex */
 	/* Go through the vertex hash to set every vertex value */
 	for (i = 0; i < vc->hash_size; i++) {
 		for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+			double *col = NULL;
 
 			if (vx->status != vtx_norm && vx->addvtx)
 				error ("Found vertex that is both deleted and cause of added bxcell");
 
 			if (doadded && vx->addvtx)				/* Cause of added bxcell */
-				vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, green);
+				col = green;
 			else if (dopres && vx->pres)			/* Preserved vertex */
-				vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, yellow);
+				col = yellow;
 			else if (dodeleted && (vx->status == vtx_sha || vx->status == vtx_del))
-				vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, red);
+				col = red;
 			else if (dooil && vx->status == vtx_oil)
-				vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, blue);
+				col = blue;
 			else if (vx->status == vtx_norm)
-				vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, white);
+				col = white;
+
+			if (col != NULL) {
+				vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->VV, col);
+			}
 		}
 	}
 
@@ -12304,26 +12935,26 @@ assdire *edgdir		/* Edge lookup for vertex */
 					col = white;
 
 				for (f = 0; f < fdi; f++)
-					vv[f] = vx->v[f] + off;
+					vv[f] = vx->VV[f] + off;
 				lix[0] = wrl->add_vertex(wrl, 2, vv);
 				for (f = 0; f < fdi; f++)
-					vv[f] = vx->v[f] - off;
+					vv[f] = vx->VV[f] - off;
 				lix[1] = wrl->add_vertex(wrl, 2, vv);
 				wrl->add_col_line(wrl, 2, lix, col);
 
 				for (f = 0; f < fdi; f++)
-					vv[f] = vx->v[f] + ((f & 1) ? off : -off);
+					vv[f] = vx->VV[f] + ((f & 1) ? off : -off);
 				lix[0] = wrl->add_vertex(wrl, 2, vv);
 				for (f = 0; f < fdi; f++)
-					vv[f] = vx->v[f] - ((f & 1) ? off : -off);
+					vv[f] = vx->VV[f] - ((f & 1) ? off : -off);
 				lix[1] = wrl->add_vertex(wrl, 2, vv);
 				wrl->add_col_line(wrl, 2, lix, col);
 
 				for (f = 0; f < fdi; f++)
-					vv[f] = vx->v[f] + ((f & 2) ? off : -off);
+					vv[f] = vx->VV[f] + ((f & 2) ? off : -off);
 				lix[0] = wrl->add_vertex(wrl, 2, vv);
 				for (f = 0; f < fdi; f++)
-					vv[f] = vx->v[f] - ((f & 2) ? off : -off);
+					vv[f] = vx->VV[f] - ((f & 2) ? off : -off);
 				lix[1] = wrl->add_vertex(wrl, 2, vv);
 				wrl->add_col_line(wrl, 2, lix, col);
 			}
@@ -12396,7 +13027,7 @@ int bxix			/* Index of bx cell causing touches */
 	double white[3] = { 1.0, 1.0, 1.0 };
 	double red[3] = { 0.8, 0.1, 0.1 };
 
-	wrl = new_vrml("add_touch_bxcells", 0, vrml_lab);
+	wrl = new_vrml("add_touch_bxcells", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
 
 	/* Gamut center point marker */
 	wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
@@ -12461,7 +13092,7 @@ int dowait				/* Wait for a return key */
 	double grey[3] = { 0.5, 0.5, 0.5 };
 	double white[3] = { 1.0, 1.0, 1.0 };
 
-	wrl = new_vrml("thinned_fwcells", 0, vrml_lab);
+	wrl = new_vrml("thinned_fwcells", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
 	wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
 
 	if (dofclabels) {