state from 2014-09-07 13:30 MESZ

9 files changed, 752 insertions, 215 deletions

diff --git a/rspl/Jamfile b/rspl/Jamfile
index a3f2e8c..dbca6d5 100644
--- a/rspl/Jamfile
+++ b/rspl/Jamfile
@@ -3,8 +3,8 @@
 
 # Optimization and Debug flags
 
-PREF_CCFLAGS 	+= $(CCOPTFLAG) ;		# Turn optimisation on
-#PREF_CCFLAGS	+= $(CCDEBUGFLAG) ;		# Debugging flags
+#PREF_CCFLAGS 	+= $(CCOPTFLAG) ;		# Turn optimisation on
+PREF_CCFLAGS	+= $(CCDEBUGFLAG) ;		# Debugging flags
 #PREF_CCFLAGS 	+= $(CCPROFFLAG) ;		# Profile flags
 #PREF_LINKFLAGS	+= $(LINKPROFFLAG) ;    # Profile flags
 #PREF_CCFLAGS		+= $(CCHEAPDEBUG) ;		# Heap Debugging flags
@@ -27,9 +27,9 @@ HDRS = ../h ../numlib ../plot $(TIFFINC) ;
 LINKLIBS = librspl ../numlib/libnum ../plot/libplot ../plot/libvrml ../icc/libicc $(TIFFLIB) $(JPEGLIB) ;
 
 # Test programs
-MainsFromSources revbench.c c1.c c1df.c t2d.c t2ddf.c t3d.c t3ddf.c tnd.c trnd.c ;
+MainsFromSources revbench.c c1.c cw1.c cw3.c c1df.c t2d.c t2ddf.c t3d.c t3ddf.c tnd.c trnd.c ;
 
-BUILD_TESTS = true ;
+BUILD_TESTS = false ;
 
 if $(BUILD_TESTS) {
 
diff --git a/rspl/Makefile.am b/rspl/Makefile.am
deleted file mode 100644
index cfa0ec2..0000000
--- a/rspl/Makefile.am
+++ /dev/null
@@ -1,15 +0,0 @@
-include $(top_srcdir)/Makefile.shared
-
-privatelib_LTLIBRARIES = librspl.la
-privatelibdir = $(pkglibdir)
-
-librspl_la_SOURCES = rspl.h rspl_imp.h mlbs.h rspl.c scat.c rev.c	\
-	rev.h gam.c spline.c opt.c
-librspl_la_LIBADD = ../numlib/libargyllnum.la ../plot/libvrml.la
-
-LDADD = ./librspl.la ../numlib/libargyllnum.la ../plot/libplot.la	\
-	../plot/libvrml.la $(X_LIBS) $(TIFF_LIBS) $(ICC_LIBS)
-
-check_PROGRAMS = revbench c1 c1df t2d t2ddf t3d t3ddf tnd trnd
-
-EXTRA_DIST = License.txt Readme.txt
diff --git a/rspl/afiles b/rspl/afiles
index af20c44..7d9c29f 100644
--- a/rspl/afiles
+++ b/rspl/afiles
@@ -16,6 +16,8 @@ gam.h
 gam.c
 revbench.c
 c1.c
+cw1.c
+cw3.c
 c1df.c
 spline.c
 t2d.c
diff --git a/rspl/c1.c b/rspl/c1.c
index 8172cef..2717920 100644
--- a/rspl/c1.c
+++ b/rspl/c1.c
@@ -19,7 +19,7 @@
 #undef DIAG
 #undef DIAG2
 #undef GLOB_CHECK
-#undef RES2			/* Do multiple test at various resolutions */
+#define RES2			/* Do multiple test at various resolutions */
 #define AVGDEV 0.0		/* Average deviation of function data */
 
 #include <stdio.h>
@@ -285,7 +285,6 @@ int main(int argc, char *argv[]) {
 				rss->fit_rspl(rss,
 			           0 | (twopass ? RSPL_2PASSSMTH : 0)
 			             | (extra   ? RSPL_EXTRAFIT2 : 0) ,
-			           0,
 			           test_points,			/* Test points */
 			           pnts,	/* Number of test points */
 			           low, high, gres,		/* Low, high, resolution of grid */
diff --git a/rspl/cw1.c b/rspl/cw1.c
new file mode 100644
index 0000000..75a4f55
--- /dev/null
+++ b/rspl/cw1.c
@@ -0,0 +1,308 @@
+
+/************************************************/
+/* Investigate various curve approximations     */
+/************************************************/
+
+/* Discrete regularized spline versions */
+/* Test variable grid spacing in 1D
+
+/* Author: Graeme Gill
+ * Date:   4/10/95
+ * Date:   5/4/96
+ *
+ * Copyright 1995, 2013 Graeme W. Gill
+ *
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ */
+
+#undef DIAG
+#undef DIAG2
+#undef GLOB_CHECK
+#define AVGDEV 0.0		/* Average deviation of function data */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <fcntl.h>
+#include <math.h>
+#include "copyright.h"
+#include "aconfig.h"
+#include "numlib.h"
+#include "plot.h"
+#include "rspl.h"
+
+double lin();
+void usage(void);
+
+#define TRIALS 1	/* Number of random trials */
+#define XRES 100		/* Plotting res */
+
+
+#define PNTS1 10
+#define GRES1 100
+
+/* Function that we're approximating */
+static double func(double val) {
+	double out;
+	int sgn = 0;
+
+	if (val < 0.0) {
+		sgn = 1;
+		val = -val;
+	}
+
+	out = pow(val, 2.5);
+
+	if (sgn)
+		out = -out;
+
+	return out;
+}
+
+/* Grid spacing power */
+#define GRIDPOW 1.5
+
+/* Scattered data sample spacing power */
+#define SAMPPOW 1.3
+
+
+co test_points[PNTS1];
+
+double ipos[GRES1];
+double *iposes[1] = { ipos };
+
+double powlike(double vv, double pp);
+
+void usage(void) {
+	fprintf(stderr,"Test 1D rspl interpolation variable grid spacing\n");
+	fprintf(stderr,"Author: Graeme W. Gill\n");
+	fprintf(stderr,"usage: c1 [options]\n");
+	fprintf(stderr," -s smooth     Use given smoothness (default 1.0)\n");
+	exit(1);
+}
+
+int main(int argc, char *argv[]) {
+	int fa,nfa;				/* argument we're looking at */
+	int i,j, n;
+	double x;
+	double xx1[XRES];
+	double yy1[6][XRES];
+	double xx2[XRES];
+	double yy2[6][XRES];
+	rspl *rss;		/* incremental solution version */
+	datai low,high;
+	int gres[MXDI];
+	double smooth = 1.0;
+	double avgdev[MXDO];
+
+	low[0] = 0.0;
+	high[0] = 1.0;
+	avgdev[0] = AVGDEV;
+
+	error_program = "c1";
+	check_if_not_interactive();
+
+	/* 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();
+
+			/* smoothness */
+			else if (argv[fa][1] == 's' || argv[fa][1] == 'S') {
+				fa = nfa;
+				if (na == NULL) usage();
+				smooth = atof(na);
+			}
+			else 
+				usage();
+		} else
+			break;
+	}
+
+	/* Only one trial */
+	for (n = 0; n < TRIALS; n++) {
+		double lrand = 0.0;		/* Amount of level randomness */
+		int pnts;
+		int fres;
+
+		pnts = PNTS1;
+		fres = GRES1; 
+		gres[0] = fres;
+
+		/* Create the object */
+		rss =  new_rspl(RSPL_NOFLAGS,
+		                1,				/* di */
+		                1);				/* fdi */
+
+		/* For each grid spacing test */
+		for (j = 0; j < 2; j++) {
+
+			for (i = 0; i < pnts; i++) {
+				double val, out;
+
+				val = i/(pnts - 1.0);
+
+				if (j == 0) {			/* Linear */
+//					out = func(val + d_rand(-0.001, 0.001));
+					out = func(val);
+					test_points[i].p[0] = val;	
+					test_points[i].v[0] = out;
+
+				} else {					/* power 2.0 grid spacing */
+					val = pow(val, SAMPPOW);	/* sample point spacing */ 
+
+//					out = func(val + d_rand(-0.001, 0.001));
+					out = func(val);
+					val = powlike(val, 1.0/GRIDPOW);	/* To grid spacing */ 
+					test_points[i].p[0] = val;
+					test_points[i].v[0] = out;
+				}
+			}
+
+			/* Set grid position info */
+			for (i = 0; i < gres[0]; i++) {
+				double val, out;
+
+				val = i/(gres[0] - 1.0);
+				if (j == 1)
+					val = powlike(val, 1.0/GRIDPOW);	/* To grid spacing */ 
+				ipos[i] = val; 
+			}
+			
+			rss->fit_rspl(rss,
+		           0,
+		           test_points,			/* Test points */
+		           pnts,	/* Number of test points */
+		           low, high, gres,		/* Low, high, resolution of grid */
+		           NULL, NULL,			/* Default data scale */
+		           smooth,				/* Smoothing */
+		           avgdev,				/* Average deviation */
+		           j == 0 ? NULL : iposes);				/* iwidth */
+
+			/* Show the result full scale */
+			for (i = 0; i < XRES; i++) {
+				co tp;	/* Test point */
+				x = i/(double)(XRES-1);
+				xx1[i] = x;
+				yy1[0][i] = func(x);
+				tp.p[0] = x;
+				if (j == 1)
+					tp.p[0] = powlike(tp.p[0], 1.0/GRIDPOW);	/* Grid spacing input */ 
+				rss->interp(rss, &tp);
+				yy1[1+j][i] = tp.v[0];
+				if (yy1[1+j][i] < -0.2)
+					yy1[1+j][i] = -0.2;
+				else if (yy1[1+j][i] > 1.2)
+					yy1[1+j][i] = 1.2;
+			}
+
+			/* Show the result magnified */
+			for (i = 0; i < XRES; i++) {
+				co tp;	/* Test point */
+				x = i/(double)(XRES-1);
+				x *= 0.1;
+				xx2[i] = x;
+				yy2[0][i] = func(x);
+				tp.p[0] = x;
+				if (j == 1)
+					tp.p[0] = powlike(tp.p[0], 1.0/GRIDPOW);	/* Grid spacing input */ 
+				rss->interp(rss, &tp);
+				yy2[1+j][i] = tp.v[0];
+				if (yy2[1+j][i] < -0.2)
+					yy2[1+j][i] = -0.2;
+				else if (yy2[1+j][i] > 1.2)
+					yy2[1+j][i] = 1.2;
+			}
+		}
+	
+		printf("Full scale: Black = lin, Red = even, Green = pow2\n");
+		do_plot6(xx1,yy1[0],yy1[1],yy1[2],NULL,NULL,NULL,XRES);
+
+		printf("Magnified: Black = lin, Red = even, Green = pow2\n");
+		do_plot6(xx2,yy2[0],yy2[1],yy2[2],NULL,NULL,NULL,XRES);
+	}
+	return 0;
+}
+
+
+/* A power-like function, based on Graphics Gems adjustment curve. */
+/* Avoids "toe" problem of pure power. */
+/* Adjusted so that "power" 2 and 0.5 agree with real power at 0.5 */
+
+double powlike(double vv, double pp) {
+	double tt, g;
+
+	if (pp >= 1.0) {
+		g = 2.0 * (pp - 1.0);
+		vv = vv/(g - g * vv + 1.0);
+	} else {
+		g = 2.0 - 2.0/pp;
+		vv = (vv - g * vv)/(1.0 - g * vv);
+	}
+
+	return vv;
+}
+
+
+/******************************************************************/
+/* Error/debug output routines */
+/******************************************************************/
+
+/* Next u function done with optimization */
+
+/* Structure to hold data for optimization function */
+struct _edatas {
+	rspl *rss;
+	int j;
+	}; typedef struct _edatas edatas;
+
+#ifdef GLOB_CHECK
+/* Overall Global optimization method */
+/* Definition of the optimization function handed to powell() */
+double efunc2(void *edata, double p[])
+	{
+	int j;
+	double rv;
+	rspl *rss = (rspl *)edata;
+	for (j = 0; j < rss->nig; j++)	/* Ugg */
+		rss->u[j].v = p[j];
+	rv = rss->efactor(rss);
+#ifdef DIAG2
+	/* printf("%c%e",cr_char,rv); */
+	printf("%e\n",rv);
+#endif
+	return rv;
+	}
+
+solveu(rss)
+rspl *rss;
+	{
+	int j;
+	double *cp;
+	double *s;
+
+	cp = dvector(0,rss->nig);
+	s = dvector(0,rss->nig);
+	for (j = 0; j < rss->nig; j++)	/* Ugg */
+		{
+		cp[j] = rss->u[j].v;
+		s[j] = 0.1;
+		}
+	powell(rss->nig,cp,s,1e-7,1000,efunc2,(void *)rss);
+	}
+#endif /* GLOB_CHECK */
diff --git a/rspl/cw3.c b/rspl/cw3.c
new file mode 100644
index 0000000..03c9b81
--- /dev/null
+++ b/rspl/cw3.c
@@ -0,0 +1,325 @@
+
+/************************************************/
+/* Investigate various curve approximations     */
+/************************************************/
+
+/* Discrete regularized spline versions */
+/* Test variable grid spacing in 3D
+
+/* Author: Graeme Gill
+ * Date:   28/11/2013
+ *
+ * Copyright 1995, 2013 Graeme W. Gill
+ *
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ */
+
+#undef DIAG
+#undef DIAG2
+#undef GLOB_CHECK
+#define AVGDEV 0.0		/* Average deviation of function data */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <fcntl.h>
+#include <math.h>
+#include "copyright.h"
+#include "aconfig.h"
+#include "counters.h"
+#include "numlib.h"
+#include "plot.h"
+#include "rspl.h"
+
+double lin();
+void usage(void);
+
+#define TRIALS 1	/* Number of random trials */
+#define XRES 100		/* Plotting res */
+
+
+#define PNTS1 20
+#define GRES1 33
+
+/* Function that we're approximating */
+static double func(double val) {
+	double out;
+	int sgn = 0;
+
+	if (val < 0.0) {
+		sgn = 1;
+		val = -val;
+	}
+
+	out = pow(val, 2.4);
+
+	if (sgn)
+		out = -out;
+
+	return out;
+}
+
+/* Grid spacing power */
+#define GRIDPOW 1.8
+
+/* Scattered data sample spacing power */
+#define SAMPPOW 1.3
+
+
+co test_points[PNTS1 * PNTS1 * PNTS1];
+
+double ipos[GRES1];
+double *iposes[3] = { ipos, ipos, ipos };
+
+double powlike(double vv, double pp);
+
+void usage(void) {
+	fprintf(stderr,"Test 3D rspl interpolation variable grid spacing\n");
+	fprintf(stderr,"Author: Graeme W. Gill\n");
+	fprintf(stderr,"usage: c1 [options]\n");
+	fprintf(stderr," -s smooth     Use given smoothness (default 1.0)\n");
+	exit(1);
+}
+
+int main(int argc, char *argv[]) {
+	int fa,nfa;				/* argument we're looking at */
+	int i, j, e, f, n;
+	double x;
+	double xx1[XRES];
+	double yy1[6][XRES];
+	double xx2[XRES];
+	double yy2[6][XRES];
+	rspl *rss;		/* incremental solution version */
+	datai low,high;
+	int gres[MXDI];
+	double smooth = 1.0;
+	double avgdev[MXDO];
+
+	low[0] = 0.0;
+	low[1] = 0.0;
+	low[2] = 0.0;
+	high[0] = 1.0;
+	high[1] = 1.0;
+	high[2] = 1.0;
+	avgdev[0] = AVGDEV;
+	avgdev[1] = AVGDEV;
+	avgdev[2] = AVGDEV;
+
+	error_program = "c1";
+	check_if_not_interactive();
+
+	/* 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();
+
+			/* smoothness */
+			else if (argv[fa][1] == 's' || argv[fa][1] == 'S') {
+				fa = nfa;
+				if (na == NULL) usage();
+				smooth = atof(na);
+			}
+			else 
+				usage();
+		} else
+			break;
+	}
+
+	for (n = 0; n < TRIALS; n++) {
+		double lrand = 0.0;		/* Amount of level randomness */
+		int di = 3;
+		int pnts = PNTS1;
+		int tpnts = PNTS1 * PNTS1 * PNTS1;
+		int fres;
+
+		fres = GRES1; 
+		gres[0] = fres;
+		gres[1] = fres;
+		gres[2] = fres;
+
+		/* Create the object */
+		rss =  new_rspl(RSPL_NOFLAGS,
+		                3,				/* di */
+		                3);				/* fdi */
+
+		/* For each grid spacing test */
+		for (j = 0; j < 2; j++) {
+			DCOUNT(gc, MXDIDO, di, 0, 0, pnts);
+			
+			/* Set the scattered data points (actually grid) */
+			DC_INIT(gc);
+			for (i = 0;  !DC_DONE(gc); i++) {
+
+				for (f = 0; f < 3; f++) {
+					double val, out;
+
+					val = gc[f]/(pnts - 1.0);
+
+					if (j == 0) {			/* Linear */
+//						out = func(val + d_rand(-0.001, 0.001));
+						out = func(val);
+						test_points[i].p[f] = val;	
+						test_points[i].v[f] = out;
+	
+					} else {					/* power 2.0 grid spacing */
+						val = pow(val, SAMPPOW);	/* sample point spacing */ 
+	
+//						out = func(val + d_rand(-0.001, 0.001));
+						out = func(val);
+						val = powlike(val, 1.0/GRIDPOW);	/* To grid spacing */ 
+						test_points[i].p[f] = val;
+						test_points[i].v[f] = out;
+					}
+				}
+				DC_INC(gc);
+			}
+
+			/* Set grid position info. This is the same for each input axis */
+			for (i = 0; i < gres[0]; i++) {
+				double val, out;
+
+				val = i/(gres[0] - 1.0);
+				if (j == 1)
+					val = powlike(val, 1.0/GRIDPOW);	/* To grid spacing */ 
+				ipos[i] = val; 
+			}
+			
+			printf("Fitting rspl....\n");
+			rss->fit_rspl(rss,
+		           0,
+		           test_points,			/* Test points */
+		           tpnts,				/* Number of test points */
+		           low, high, gres,		/* Low, high, resolution of grid */
+		           NULL, NULL,			/* Default data scale */
+		           smooth,				/* Smoothing */
+		           avgdev,				/* Average deviation */
+		           j == 0 ? NULL : iposes);				/* iwidth */
+
+			/* Show the result full scale of one axis and corresponding output */
+			for (i = 0; i < XRES; i++) {
+				co tp;	/* Test point */
+				x = i/(double)(XRES-1);
+				xx1[i] = x;
+				yy1[0][i] = func(x);
+				tp.p[0] = x;
+				tp.p[1] = tp.p[2] = 0.0;
+				if (j == 1)
+					tp.p[0] = powlike(tp.p[0], 1.0/GRIDPOW);	/* Grid spacing input */ 
+				rss->interp(rss, &tp);
+				yy1[1+j][i] = tp.v[0];
+				if (yy1[1+j][i] < -0.2)
+					yy1[1+j][i] = -0.2;
+				else if (yy1[1+j][i] > 1.2)
+					yy1[1+j][i] = 1.2;
+			}
+
+			/* Show the result magnified */
+			for (i = 0; i < XRES; i++) {
+				co tp;	/* Test point */
+				x = i/(double)(XRES-1);
+				x *= 0.1;
+				xx2[i] = x;
+				yy2[0][i] = func(x);
+				tp.p[0] = x;
+				tp.p[1] = tp.p[2] = 0.0;
+				if (j == 1)
+					tp.p[0] = powlike(tp.p[0], 1.0/GRIDPOW);	/* Grid spacing input */ 
+				rss->interp(rss, &tp);
+				yy2[1+j][i] = tp.v[0];
+				if (yy2[1+j][i] < -0.2)
+					yy2[1+j][i] = -0.2;
+				else if (yy2[1+j][i] > 1.2)
+					yy2[1+j][i] = 1.2;
+			}
+		}
+	
+		printf("Full scale: Black = lin, Red = even, Green = pow2\n");
+		do_plot6(xx1,yy1[0],yy1[1],yy1[2],NULL,NULL,NULL,XRES);
+
+		printf("Magnified: Black = lin, Red = even, Green = pow2\n");
+		do_plot6(xx2,yy2[0],yy2[1],yy2[2],NULL,NULL,NULL,XRES);
+	}
+	return 0;
+}
+
+/* A power-like function, based on Graphics Gems adjustment curve. */
+/* Avoids "toe" problem of pure power. */
+/* Adjusted so that "power" 2 and 0.5 agree with real power at 0.5 */
+
+double powlike(double vv, double pp) {
+	double tt, g;
+
+	if (pp >= 1.0) {
+		g = 2.0 * (pp - 1.0);
+		vv = vv/(g - g * vv + 1.0);
+	} else {
+		g = 2.0 - 2.0/pp;
+		vv = (vv - g * vv)/(1.0 - g * vv);
+	}
+
+	return vv;
+}
+
+
+/******************************************************************/
+/* Error/debug output routines */
+/******************************************************************/
+
+/* Next u function done with optimization */
+
+/* Structure to hold data for optimization function */
+struct _edatas {
+	rspl *rss;
+	int j;
+	}; typedef struct _edatas edatas;
+
+#ifdef GLOB_CHECK
+/* Overall Global optimization method */
+/* Definition of the optimization function handed to powell() */
+double efunc2(void *edata, double p[])
+	{
+	int j;
+	double rv;
+	rspl *rss = (rspl *)edata;
+	for (j = 0; j < rss->nig; j++)	/* Ugg */
+		rss->u[j].v = p[j];
+	rv = rss->efactor(rss);
+#ifdef DIAG2
+	/* printf("%c%e",cr_char,rv); */
+	printf("%e\n",rv);
+#endif
+	return rv;
+	}
+
+solveu(rss)
+rspl *rss;
+	{
+	int j;
+	double *cp;
+	double *s;
+
+	cp = dvector(0,rss->nig);
+	s = dvector(0,rss->nig);
+	for (j = 0; j < rss->nig; j++)	/* Ugg */
+		{
+		cp[j] = rss->u[j].v;
+		s[j] = 0.1;
+		}
+	powell(rss->nig,cp,s,1e-7,1000,efunc2,(void *)rss);
+	}
+#endif /* GLOB_CHECK */
diff --git a/rspl/rev.c b/rspl/rev.c
index 237dbe3..1dbf6cc 100644
--- a/rspl/rev.c
+++ b/rspl/rev.c
@@ -22,11 +22,24 @@
     Should fix the clipping case so that a direction weighting
 	funtion can be applied. This should be used just like
 	the perceptual case to increase L* constance for dark
-    colors. Will this stuff up the geometric consistency though ?
-	[ See fill_nncell(), fill_nncell() and users of calc_fwd_nn_cell_list(),
-	  ie. nnearest_clip_solve(), clipn_setsort() etc. ]
+    colors. This would entail large scale changes though,
+	since a lot of code assumes minimal euclidean distance
+	goal, from the cell selection structure [ See fill_nncell(),
+	fill_nncell() and users of calc_fwd_nn_cell_list() ] and
+	the within cell computation [ ie. See  nnearest_clip_solve(),
+	clipn_setsort() etc. ]
+	XYZ PCS couldn't work with a simple weighting - it would have
+	to be a position dependent weighting.
 	The SVD least squares computation case makes this hard to change ?
 	Would have to feed in a weighting function, or can it be general ?
+	-
+	Can this be solved some other way, ie. by using gamut
+	mapping type look up ? Problem is precision.
+	-
+	Vector clip could be used (if intent can be turned
+	into computable vector clip direction), but it is slow,
+	because it search all cells from source until it
+	hits surface.
 
 	Allow function callback to set auxiliary values for 
 	flag RSPL_AUXLOCUS. 
@@ -1437,7 +1450,8 @@ unsigned int tcount		/* grid touch count for this operation */
 							nunlk++;
 					}
 					fprintf(stdout,"Diagnostic: rev.sz = %lu, rev.max_sz = %lu, numlocked = %d, nunlk = %d\n",
-					               rc->s->rev.sz, rc->s->rev.max_sz, rc->nunlocked,nunlk);
+					               (unsigned long)rc->s->rev.sz, (unsigned long)rc->s->rev.max_sz,
+					               rc->nunlocked, nunlk);
 					error("Not enough memory to process in chunks");
 				}
 				break;		/* cache has run out of room, so abandon, and do it next time */
@@ -4971,7 +4985,7 @@ rspl *s		/* Pointer to rspl grid */
 								g_no_rev_cache_instances > 1 ? "are" : "is",
 			                    g_no_rev_cache_instances,
 								g_no_rev_cache_instances > 1 ? "s" : "",
-			                    ram_portion/1000000);
+			                    (unsigned long)ram_portion/1000000);
 		}
 	}
 
@@ -6165,7 +6179,7 @@ if (prop != NULL) {
 								g_no_rev_cache_instances > 1 ? "are" : "is",
 			                    g_no_rev_cache_instances,
 								g_no_rev_cache_instances > 1 ? "s" : "",
-			                    ram_portion/1000000);
+			                    (unsigned long)ram_portion/1000000);
 	}
 	s->rev.rev_valid = 1;
 
@@ -6233,7 +6247,7 @@ rspl *s		/* Pointer to rspl grid */
 								g_no_rev_cache_instances > 1 ? "are" : "is",
 			                    g_no_rev_cache_instances,
 								g_no_rev_cache_instances > 1 ? "s" : "",
-			                    ram_portion/1000000);
+			                    (unsigned long)ram_portion/1000000);
 		}
 	}
 	s->rev.rev_valid = 0;
@@ -6497,7 +6511,7 @@ rspl *s
 			if (g_avail_ram > safe_max_vmem) {
 				g_avail_ram = safe_max_vmem;
 				if (s->verbose && repsr == 0)
-					fprintf(stdout,"%cTrimmed maximum cache RAM to %lu Mbytes to allow for VM limit\n",cr_char,g_avail_ram/1000000);
+					fprintf(stdout,"%cTrimmed maximum cache RAM to %lu Mbytes to allow for VM limit\n",cr_char,(unsigned long)g_avail_ram/1000000);
 			}
 		}
 	
@@ -6516,7 +6530,7 @@ rspl *s
 		}
 		if (max_vmem != 0 && g_avail_ram > max_vmem && repsr == 0) {
 			g_avail_ram = (size_t)(0.95 * max_vmem);
-			fprintf(stdout,"%cARGYLL_REV_CACHE_MULT * RAM trimmed to %lu Mbytes to allow for VM limit\n",cr_char,g_avail_ram/1000000);
+			fprintf(stdout,"%cARGYLL_REV_CACHE_MULT * RAM trimmed to %lu Mbytes to allow for VM limit\n",cr_char,(unsigned long)g_avail_ram/1000000);
 		}
 	}
 
@@ -6525,7 +6539,7 @@ rspl *s
 
 	DBG(("reverse cache max memory = %d Mbytes\n",s->rev.max_sz/1000000));
 	if (s->verbose && repsr == 0) {
-		fprintf(stdout, "%cRev cache RAM = %lu Mbytes\n",cr_char,g_avail_ram/1000000);
+		fprintf(stdout, "%cRev cache RAM = %lu Mbytes\n",cr_char,(unsigned long)g_avail_ram/1000000);
 		repsr = 1;
 	}
 
diff --git a/rspl/rspl.c b/rspl/rspl.c
index 13a1776..23ef6b4 100644
--- a/rspl/rspl.c
+++ b/rspl/rspl.c
@@ -866,7 +866,7 @@ rspl *s
 
 /* ============================================ */
 /* Initialize the grid from a provided function. By default the grid */
-/* values are set to exactly the value returned fy func(), unless the */
+/* values are set to exactly the value returned by func(), unless the */
 /* RSPL_SET_APXLS flag is set, in which case an attempt is made to have */
 /* the grid points represent a least squares aproximation to the underlying */
 /* surface, by using extra samples in the middle of grid cells. */
diff --git a/rspl/scat.c b/rspl/scat.c
index b4ed978..ec7e99c 100644
--- a/rspl/scat.c
+++ b/rspl/scat.c
@@ -128,10 +128,6 @@
 /* algorithm parameters [Release defaults] */
 #undef POINTWEIGHT		/* [Undef] Increas smoothness weighting proportional to number of points */
 #define INCURVEADJ		/* [Defined] Adjust smoothness criteria for input curve grid spacing */
-#define EXTRA_SURFACE_SMOOTHING /* [Defined] Stiffen surface points to comp. for single ended. */
-						/* The following are available, but the smoothing table is */
-						/* not setup for them, and they are not sufficiently different */
-						/* from the default smoothing to be useful. */
 #define ENABLE_2PASSSMTH	/* [Define] Enable 2 pass smooth using Gaussian filter */
 #define ENABLE_EXTRAFIT	/* [Undef] Enable the extra fit option. Good to combat high smoothness. */
 
@@ -142,12 +138,12 @@
 
 /* Experimental set: */
 
-#pragma message("!!!!!!!!! Experimental config set !!!!!!!!!")
+#pragma message("!!!!!!!!! Experimental hi-accuracy config set !!!!!!!!!")
 
 #define TOL 1e-12		/* Tollerance of result - usually 1e-5 is best. */
 #define TOL_IMP 1.0		/* Minimum error improvement to continue - reduces accuracy (1.0 == off) */
 #undef GRADUATED_TOL	/* Speedup attemp - use reduced tollerance for prior grids. */
-#define GRATIO 2.0		/* Multi-grid resolution ratio */
+#define GRATIO 1.5		/* Multi-grid resolution ratio */
 #undef OVERRLX 		/* Use over relaxation factor when progress slows (worse accuracy ?) */
 #define JITTERS 0		/* Number of 1D conjugate solve itters */
 #define CONJ_TOL 1.0	/* Extra tolereance on 1D conjugate solution times TOL. */
@@ -160,8 +156,8 @@
 
 /* Release set: */
 
-#define TOL 1e-6		/* [1e-6] Tollerance of result - usually 1e-5 is best. */
-#define TOL_IMP 0.998	/* [0.998] Minimum error improvement to continue - reduces accuracy (1.0 == off) */
+#define TOL 1e-7		/* [1e-6] Tollerance of result - usually 1e-5 is best. */
+#define TOL_IMP 0.999	/* [0.998] Minimum error improvement to continue - reduces accuracy (1.0 == off) */
 #undef GRADUATED_TOL	/* [Undef] Speedup attemp - use reduced tollerance for prior grids. */
 #define GRATIO 2.0		/* [2.0] Multi-grid resolution ratio */
 #undef OVERRLX 			/* [Undef] Use over relaxation when progress slows (worse accuracy ?) */
@@ -262,7 +258,7 @@ static void free_cj_arrays(cj_arrays *ta);
 static int add_rspl_imp(rspl *s, int flags, void *d, int dtp, int dno);
 static mgtmp *new_mgtmp(rspl *s, int gres[MXDI], int f);
 static void free_mgtmp(mgtmp *m);
-static void setup_solve(mgtmp *m, int final);
+static void setup_solve(mgtmp *m);
 static void solve_gres(mgtmp *m, cj_arrays *ta, double tol, int final);
 static void init_soln(mgtmp  *m1, mgtmp  *m2);
 static void comp_ccv(mgtmp *m);
@@ -656,8 +652,7 @@ add_rspl_imp(
 					if (s->tpsm && s->tpsm2 != 0) {	/* 2nd pass of 2 pass smoothing */
 						init_ccv(m);				/* Downsample m->ccv from s->g.ccv */
 					}
-//					setup_solve(m, nn == (s->niters-1));
-					setup_solve(m, 1);
+					setup_solve(m);
 
 					if (nn == 0) {						/* Make sure we have an initial x[] */
 						for (i = 0; i <  m->g.no; i++)
@@ -921,9 +916,20 @@ free_data(rspl *s) {
 /* to them, would have an average sample deviation of 0.05. */
 /* For normally distributed errors, the average deviation is */
 /* aproximately 0.564 times the standard deviation. (0.564 * sqrt(variance)) */
-/* This table is appropriate for the default rspl algorithm + slight EXTRA_SURFACE_SMOOTHING, */
+/* This table is appropriate for the default rspl algorithm, */
 /* and is NOT setup for RSPL_2PASSSMTH or RSPL_EXTRAFIT2 !! */
 /* SMOOTH */
+
+/* There are still issues with all this - the level of smoothing actually */
+/* depends on the degree of fit of the underlying model - ie. how close */
+/* to straight the mapping is. To get actual noise reduction under these */
+/* conditions is harder than when there is some curvature to "tension" things. */
+/* This is evident is Lab vs. XYZ display profiles, and there is code */
+/* in xlut.c that tries to adjust to this. */
+
+/* !!! Possible answer - should be using third order differences */
+/*     for controlling smoothness, not second order (curvature) ?? */
+
 // ~~99
 static double opt_smooth(
 	rspl *s,
@@ -963,7 +969,7 @@ static double opt_smooth(
 	
 
 	/* New for V1.10, from smtmpp using sRGB, EpsonR1800, Hitachi2112, */
-	/* Fogra39L, Canon1180, Epson10K, with low EXTRA_SURFACE_SMOOTHING. */
+	/* Fogra39L, Canon1180, Epson10K (did use EXTRA_SURFACE_SMOOTHING). */
 
 	/* Main lookup table, by [di][ncix][adix]: */
 	/* Values are log of smoothness value. */
@@ -1349,8 +1355,7 @@ static void free_mgtmp(mgtmp  *m) {
  */
 
 static void setup_solve(
-	mgtmp  *m,		/* initialized grid temp structure */
-	int final		/* nz if final resolution (activate EXTRA_SURFACE_SMOOTHING) */
+	mgtmp  *m		/* initialized grid temp structure */
 ) {
 	rspl *s = m->s;
 	int di   = s->di;
@@ -1470,7 +1475,6 @@ static void setup_solve(
 			b[i] = 0.0;
 	}
 
-#ifdef NEVER
 	/* Production version, without extra edge weight */
 
 	/* Accumulate curvature dependent factors to the triangular A matrix. */
@@ -1499,20 +1503,46 @@ static void setup_solve(
 			ki-1 *   w[i-1]         *   w[i-1] * u[i-1]
 			ki   * -(w[i-1] + w[i]) * -(w[i-1] + w[i]) * u[i])
 			ki+1 *   w[i]           *   w[i] * u[i+1]
+
+		The weighting is adjusted to compensate for increased
+		grid density by the dw weight, and for the different
+		scaling either side of a grid point with the w0 and w1 scale.
+		Both are normalised so that the overall smoothness scaling
+		factor applies.
 	 */
 
-	{		/* Setting this up from scratch */
+	{	/* Setting this up from scratch */
+		double dwtw[MXDIDO];		/* Density weight normalizer */
+		int wtwc[MXDIDO];
 		ECOUNT(gc, MXDIDO, di, 0, gres, 0);
-		EC_INIT(gc);
 
 		for (oawt = 0.0, i = 1; i < 21; i++)
 			oawt += wvals[i];
 		oawt *= wvals[0];
 
+#define EN_DW			/* Enable grid density weighting */
+#define EN_CW			/* Enable grid scale weighting */
+
+		/* Compute the ipos[] weight normalisation factors */
+		for (e = 0; e < di; e++) {
+			if (m->g.ipos[e] == NULL)
+				continue;
+			dwtw[e] = 0.0;
+			for (i = 1; i < gres[e]; i++) {
+				double w;
+				w = fabs(m->g.ipos[e][i] - m->g.ipos[e][i-1]);
+//printf("[%d][%d] w = %f\n",e,i,w);
+				dwtw[e] += w;
+			}
+			dwtw[e] /= (gres[e] - 1.0);		/* Average weights */
+//printf("dwtw[%d] = %f\n",e,dwtw[e]);
+		}
+
+		EC_INIT(gc);
 		for (i = 0; i < gno; i++) {
 
 			for (e = 0; e < di; e++) {	/* For each curvature direction */
-				double w0, w1, tt;
+				double dw, w0, w1, tt;
 				double cw = 2.0 * m->sf.cw[e];	/* Overall curvature weight */
 				cw *= s->d.vw[f];				/* Scale curvature weight for data range */
 
@@ -1521,14 +1551,24 @@ static void setup_solve(
 				if ((gc[e]-2) >= 0) {
 					/* double kw = cw * gp[UO_C(e,1)].k; */	/* Cell bellow k value */
 					double kw = cw;
-					w0 = w1 = 1.0;
+					dw = w1 = 1.0;
 					if (m->g.ipos[e] != NULL) {
+
 						w0 = fabs(m->g.ipos[e][gc[e]-1] - m->g.ipos[e][gc[e]-2]);
 						w1 = fabs(m->g.ipos[e][gc[e]-0] - m->g.ipos[e][gc[e]-1]);
-						tt = sqrt(w0 * w1);		/* Normalise overall width weighting effect */
+//printf("raw [%d][%d] w0 = %f, w1 = %f\n",gc[e],i,w0,w1);
+#ifdef EN_DW
+						dw = 0.5 * (w0 + w1);			/* in this direction around -1  */
+						dw = dw/dwtw[e];				/* normalise */
+#endif
+						tt = sqrt(w0 * w1);
 						w1 = tt/w1;
+#ifndef EN_CW
+						w1 = 1.0;
+#endif
+//printf("[%d][%d] dw = %f, w1 = %f\n",gc[e],i,dw,w1);
 					}
-					A[i][ixcol[0]] += w1 * w1 * kw;
+					A[i][ixcol[0]] += dw * w1 * w1 * kw;
 					if (ccv != NULL)
 						b[i] += kw * (w1) * ccv[i - gci[e]][e]; /* Curvature compensation value */
 				}
@@ -1537,16 +1577,25 @@ static void setup_solve(
 				if ((gc[e]-1) >= 0 && (gc[e]+1) < gres[e]) {
 					/* double kw = cw * gp->k;  */		/* This cells k value */
 					double kw = cw; 
-					w0 = w1 = 1.0;
+					dw = w0 = w1 = 1.0;
 					if (m->g.ipos[e] != NULL) {
 						w0 = fabs(m->g.ipos[e][gc[e]-0] - m->g.ipos[e][gc[e]-1]);
 						w1 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]-0]);
+//printf("raw [%d][%d] w0 = %f, w1 = %f\n",gc[e],i,w0,w1);
+#ifdef EN_DW
+						dw = 0.5 * (w0 + w1);			/* in this direction around 0  */
+						dw = dw/dwtw[e];				/* normalise */
+#endif
 						tt = sqrt(w0 * w1);
 						w0 = tt/w0;
 						w1 = tt/w1;
+#ifndef EN_CW
+						w0 = w1 = 1.0;
+#endif
+//printf("[%d][%d] dw = %f, w0 = %f, w1 = %f\n",gc[e],i,dw,w0,w1);
 					}
-					A[i][ixcol[0]]      += -(w0 + w1) * -(w0 + w1) * kw;
-					A[i][ixcol[gci[e]]] += -(w0 + w1) * w1 * kw * oawt;
+					A[i][ixcol[0]]      += dw * -(w0 + w1) * -(w0 + w1) * kw;
+					A[i][ixcol[gci[e]]] += dw * -(w0 + w1) * w1 * kw * oawt;
 					if (ccv != NULL)
 						b[i] += kw * -(w0 + w1) * ccv[i][e]; /* Curvature compensation value */
 				}
@@ -1555,173 +1604,33 @@ static void setup_solve(
 				if ((gc[e]+2) < gres[e]) {
 					/* double kw = cw * gp[UO_C(e,2)].k;	*/ /* Cell above k value */
 					double kw = cw;
-					w0 = w1 = 1.0;
+					dw = w0 = w1 = 1.0;
 					if (m->g.ipos[e] != NULL) {
 						w0 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]+0]);
 						w1 = fabs(m->g.ipos[e][gc[e]+2] - m->g.ipos[e][gc[e]+1]);
-						tt = sqrt(w0 * w1);
-						w0 = tt/w0;
-						w1 = tt/w1;
-					}
-					A[i][ixcol[0]]          += w0 * w0 * kw;
-					A[i][ixcol[gci[e]]]     += w0 * -(w0 + w1) * kw;
-					A[i][ixcol[2 * gci[e]]] += w0 * w1 * kw;
-					if (ccv != NULL)
-						b[i] += kw * -(w0 + w1) * ccv[i][e]; /* Curvature compensation value */
-				}
-			}
-			EC_INC(gc);
-		}
-	}
-#endif /* NEVER */
-
-#ifdef ALWAYS
-	/* Production version that allows for extra weight on grid edges */
-
-	/* Accumulate curvature dependent factors to the triangular A matrix. */
-	/* Because it's triangular, we compute and add in all the weighting */
-	/* factors at and to the right of each cell. */
-
-	/* The ipos[] factor is to allow for the possibility that the */
-	/* grid spacing may be non-uniform in the colorspace where the */
-	/* function being modelled is smooth. Our curvature computation */
-	/* needs to make allowsance for this fact in computing the */
-	/* node value differences that equate to zero curvature. */ 
-	/*
-		The old curvature fixed grid spacing equation was:
-			ki * (u[i-1] - 2 * u[i] + u[i+1])^2
-		with derivatives wrt each node:
-			ki-1 *  1 * 2 * u[i-1] 
-			ki   * -2 * 2 * u[i]
-			ki+1 *  1 * 2 * u[i+1]
-
-		Allowing for scaling of each grid difference by w[i-1] and w[i],
-		where w[i-1] corresponds to the width of cell i-1 to i,
-		and w[i] corresponds to the width of cell i to i+1:
-			ki * (w[i-1] * (u[i-1] - u[i]) + w[i] * (u[i+1] - u[i[))^2
-		=	ki * (w[i-1] * u[i-1] - (w[i-1] + w[i]) * u[i]) + w[i] * u[i+1])^2
-		with derivatives wrt each node:
-			ki-1 *   w[i-1]         *   w[i-1] * u[i-1]
-			ki   * -(w[i-1] + w[i]) * -(w[i-1] + w[i]) * u[i])
-			ki+1 *   w[i]           *   w[i] * u[i+1]
-	 */
-	{		/* Setting this up from scratch */
-		ECOUNT(gc, MXDIDO, di, 0, gres, 0);
-#ifdef EXTRA_SURFACE_SMOOTHING
-//		double  k0w = 4.0, k1w = 1.3333;	/* Extra stiffness */
-//		double  k0w = 3.0, k1w = 1.26;		/* Some extra stiffness */
-		double  k0w = 2.0, k1w = 1.15;		/* A little extra stiffness */
-#else
-		double  k0w = 1.0, k1w = 1.0;		/* No extra weights */
+//printf("raw [%d][%d] w0 = %f, w1 = %f\n",gc[e],i,w0,w1);
+#ifdef EN_DW
+						dw = 0.5 * (w0 + w1);			/* in this direction around +1  */
+						dw = dw/dwtw[e];				/* normalise */
 #endif
-
-		EC_INIT(gc);
-		for (oawt = 0.0, i = 1; i < 21; i++)
-			oawt += wvals[i];
-		oawt *= wvals[0];
-
-		if (final == 0)
-			k0w = k1w = 1.0;			/* Activate extra edge smoothing on final grid ? */
-
-		for (i = 0; i < gno; i++) {
-			int k;
-
-			/* We're creating the equation cooeficients for solving the */
-			/* partial derivative equation w.r.t. node point i. */
-			/* Due to symetry in the smoothness interactions, only */
-			/* the triangle cooeficients of neighbour nodes is needed. */
-			for (e = 0; e < di; e++) {	/* For each curvature direction */
-				double kw, w0, w1, tt;
-				double cw = 2.0 * m->sf.cw[e];	/* Overall curvature weight */
-				double xx = 1.0;				/* Extra edge weighing */
-				cw *= s->d.vw[f];				/* Scale curvature weight for data range */
-
-				/* weight factor for outer or 2nd outer in other dimensions */
-				for (k = 0; k < di; k++) {
-					if (k == e)
-						continue;
-					if (gc[k] == 0 || gc[k] == (gres[k]-1))
-						xx *= k0w;
-					else if (gc[k] == 1 || gc[k] == (gres[k]-2))
-						xx *= k1w;
-				}
-
-				/* If at least two above lower edge in this dimension */
-				/* Add influence on Curvature of cell below */
-				if ((gc[e]-2) >= 0) {
-					/* double kw = cw * gp[-gc[e]].k; */	/* Cell bellow k value */
-					kw = cw * xx;
-					w0 = w1 = 1.0;
-					if (m->g.ipos[e] != NULL) {
-						w0 = fabs(m->g.ipos[e][gc[e]-1] - m->g.ipos[e][gc[e]-2]);
-						w1 = fabs(m->g.ipos[e][gc[e]-0] - m->g.ipos[e][gc[e]-1]);
-						tt = sqrt(w0 * w1);		/* Normalise overall width weighting effect */
-						w1 = tt/w1;
-					}
-					if ((gc[e]-2) == 0 || (gc[e]+0) == (gres[e]-1))
-						kw *= k0w;
-					else if ((gc[e]-2) == 1 || (gc[e]-0) == (gres[e]-2))
-						kw *= k1w;
-					A[i][ixcol[0]] += kw * (w1) * w1;
-					if (ccv != NULL) {
-//printf("~1 tweak b[%d] by %e\n",i,kw * (w1) * ccv[i - gci[e]][e]);
-						b[i] += kw * (w1) * ccv[i - gci[e]][e]; /* Curvature compensation value */
-					}
-				}
-				/* If not one from upper or lower edge in this dimension */
-				/* Add influence on Curvature of this cell */
-				if ((gc[e]-1) >= 0 && (gc[e]+1) < gres[e]) {
-					/* double kw = cw * gp->k;  */		/* This cells k value */
-					kw = cw * xx;
-					w0 = w1 = 1.0;
-					if (m->g.ipos[e] != NULL) {
-						w0 = fabs(m->g.ipos[e][gc[e]-0] - m->g.ipos[e][gc[e]-1]);
-						w1 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]-0]);
 						tt = sqrt(w0 * w1);
 						w0 = tt/w0;
 						w1 = tt/w1;
+#ifndef EN_CW
+						w0 = w1 = 1.0;
+#endif
+//printf("[%d][%d] dw = %f, w0 = %f, w1 = %f\n",gc[e],i,dw,w0,w1);
 					}
-					if ((gc[e]-1) == 0 || (gc[e]+1) == (gres[e]-1))
-						kw *= k0w;
-					else if ((gc[e]-1) == 1 || (gc[e]+1) == (gres[e]-2))
-						kw *= k1w;
-					A[i][ixcol[0]]      += kw * -(w0 + w1) * -(w0 + w1);
-					A[i][ixcol[gci[e]]] += kw * -(w0 + w1) * w1 * oawt;
-					if (ccv != NULL) {
-//printf("~1 tweak b[%d] by %e\n",i, kw * -(w0 + w1) * ccv[i][e]);
+					A[i][ixcol[0]]          += dw * w0 * w0 * kw;
+					A[i][ixcol[gci[e]]]     += dw * w0 * -(w0 + w1) * kw;
+					A[i][ixcol[2 * gci[e]]] += dw * w0 * w1 * kw;
+					if (ccv != NULL)
 						b[i] += kw * -(w0 + w1) * ccv[i][e]; /* Curvature compensation value */
-					}
-				}
-				/* If at least two below the upper edge in this dimension */
-				/* Add influence on Curvature of cell above */
-				if ((gc[e]+2) < gres[e]) {
-					/* double kw = cw * gp[gc[e]].k;	*/ /* Cell above k value */
-					kw = cw * xx;
-					w0 = w1 = 1.0;
-					if (m->g.ipos[e] != NULL) {
-						w0 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]+0]);
-						w1 = fabs(m->g.ipos[e][gc[e]+2] - m->g.ipos[e][gc[e]+1]);
-						tt = sqrt(w0 * w1);
-						w0 = tt/w0;
-						w1 = tt/w1;
-					}
-					if ((gc[e]+0) == 0 || (gc[e]+2) == (gres[e]-1))
-						kw *= k0w;
-					else if ((gc[e]+0) == 1 || (gc[e]+2) == (gres[e]-2))
-						kw *= k1w;
-					A[i][ixcol[0]]          += kw * (w0) * w0;
-					A[i][ixcol[gci[e]]]     += kw * (w0) * -(w0 + w1);
-					A[i][ixcol[2 * gci[e]]] += kw * (w0) * w1;
-					if (ccv != NULL) {
-//printf("~1 tweak b[%d] by %e\n",i, kw * (w0) * ccv[i + gci[e]][e]);
-						b[i] += kw * (w0) * ccv[i + gci[e]][e]; /* Curvature compensation value */
-					}
 				}
 			}
 			EC_INC(gc);
 		}
 	}
-#endif /* ALWAYS */
 
 #ifdef DEBUG
 	printf("After adding curvature equations:\n");
@@ -1736,7 +1645,6 @@ static void setup_solve(
 
 	nbsum = 0.0;	/* Zero sum of b[] squared */
 
-#ifdef ALWAYS
 	/* Accumulate weak default function factors. These are effectively a */
 	/* weak "data point" exactly at each grid point. */
 	/* (Note we're not currently doing this in a cache friendly order,   */
@@ -1775,9 +1683,7 @@ static void setup_solve(
 #endif /* DEBUG */
 
 	}
-#endif /* ALWAYS */
 
-#ifdef ALWAYS
 	/* Accumulate data point dependent factors */
 	for (n = 0; n < dno; n++) {		/* Go through all the data points */
 		int j,k;
@@ -1815,8 +1721,6 @@ static void setup_solve(
 		nbsum = 1e-4;
 	m->q.normb = nbsum;
 
-#endif /* ALWAYS */ 
-
 #ifdef DEBUG
 	printf("After adding data point equations:\n");
 	for (i = 0; i < gno; i++) {