summaryrefslogtreecommitdiff
path: root/rspl
diff options
context:
space:
mode:
authorJörg Frings-Fürst <debian@jff-webhosting.net>2015-05-01 16:24:15 +0200
committerJörg Frings-Fürst <debian@jff-webhosting.net>2015-05-01 16:24:15 +0200
commita30ba67504ffd12c4db499adbb5ce47a7d1f6036 (patch)
tree9ae1a7e3849dda6bbb5c578232f6f2fa5b2e7e7e /rspl
parent89e99e8a827859729729dfc92d74be4a8f96f1a4 (diff)
parent094535c010320967639e8e86f974d878e80baa72 (diff)
New release 1.7.0
Diffstat (limited to 'rspl')
-rw-r--r--rspl/Jamfile15
-rw-r--r--rspl/c1.c45
-rw-r--r--rspl/c1df.c3
-rw-r--r--rspl/cw1.c9
-rw-r--r--rspl/cw3.c5
-rw-r--r--rspl/gam.c35
-rw-r--r--rspl/revbench.c1
-rw-r--r--rspl/rspl.c149
-rw-r--r--rspl/rspl.h38
-rw-r--r--rspl/rspl1.h2
-rw-r--r--rspl/scat.c1597
-rw-r--r--rspl/smtmpp.c323
-rw-r--r--rspl/smtnd.c398
-rw-r--r--rspl/t2d.c21
-rw-r--r--rspl/t2ddf.c1
-rw-r--r--rspl/t3d.c21
-rw-r--r--rspl/t3ddf.c17
-rw-r--r--rspl/tnd.c5
18 files changed, 1602 insertions, 1083 deletions
diff --git a/rspl/Jamfile b/rspl/Jamfile
index dbca6d5..37c2c37 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
@@ -24,7 +24,7 @@ Headers = rspl.h ;
Library librspl : rspl.c $(SCAT).c rev.c gam.c spline.c opt.c : : : ../h ../numlib ../plot ;
HDRS = ../h ../numlib ../plot $(TIFFINC) ;
-LINKLIBS = librspl ../numlib/libnum ../plot/libplot ../plot/libvrml ../icc/libicc $(TIFFLIB) $(JPEGLIB) ;
+LINKLIBS = librspl ../plot/libplot ../numlib/libnum ../numlib/libui ../plot/libvrml ../icc/libicc $(TIFFLIB) $(JPEGLIB) ;
# Test programs
MainsFromSources revbench.c c1.c cw1.c cw3.c c1df.c t2d.c t2ddf.c t3d.c t3ddf.c tnd.c trnd.c ;
@@ -36,7 +36,7 @@ if $(BUILD_TESTS) {
HDRS = ../h ../numlib ../plot ../icc ../rspl ../xicc ../gamut ../cgats ../spectro $(TIFFINC) ;
LINKLIBS = ../xicc/libxicc ../gamut/libgamut ../spectro/libinsttypes librspl
../cgats/libcgats ../icc/libicc ../plot/libplot ../plot/libvrml
- ../numlib/libnum $(TIFFLIB) $(JPEGLIB) ;
+ ../numlib/libnum ../numlib/libui $(TIFFLIB) $(JPEGLIB) ;
# Smoothness factor tuning test in Nd.
Main smtnd : smtnd.c ;
@@ -58,8 +58,13 @@ if $(BUILD_TESTS) {
}
-# Main tt : tt.c : : ../xicc : : : ../plot/libvrml ../icc/libicc ;
+# Test code
+if $(HOME) = "d:\\usr\\graeme" && $(PWD) = "/src/argyll/rspl" {
+ # Main tt : tt.c : : ../xicc : : : ../plot/libvrml ../icc/libicc ;
+ Main crossv : crossv.c : : : ../numlib ../plot : : ../plot/libplot ../numlib/libui ../numlib/libnum ;
+
+}
if $(BUILD_JUNK) {
HDRS = ../h ../numlib ;
diff --git a/rspl/c1.c b/rspl/c1.c
index 2717920..31849ee 100644
--- a/rspl/c1.c
+++ b/rspl/c1.c
@@ -20,7 +20,7 @@
#undef DIAG2
#undef GLOB_CHECK
#define RES2 /* Do multiple test at various resolutions */
-#define AVGDEV 0.0 /* Average deviation of function data */
+#define AVGDEV 0.005 /* Average deviation of function data */
#include <stdio.h>
#include <stdlib.h>
@@ -29,10 +29,10 @@
#include "copyright.h"
#include "aconfig.h"
#include "numlib.h"
-#include "plot.h"
#include "rspl.h"
+#include "plot.h"
+#include "ui.h"
-double lin();
void usage(void);
#define TRIALS 20 /* Number of random trials */
@@ -87,8 +87,7 @@ void usage(void) {
fprintf(stderr,"Author: Graeme W. Gill\n");
fprintf(stderr,"usage: c1 [options]\n");
fprintf(stderr," -s smooth Use given smoothness (default 1.0)\n");
- fprintf(stderr," -2 Use two pass smoothing\n");
- fprintf(stderr," -x Use extra fitting\n");
+ fprintf(stderr," -x Use auto smoothing\n");
exit(1);
}
@@ -102,8 +101,7 @@ int main(int argc, char *argv[]) {
datai low,high;
int gres[MXDI];
double smooth = 1.0;
- int twopass = 0;
- int extra = 0;
+ int autosm = 0;
double avgdev[MXDO];
low[0] = 0.0;
@@ -134,17 +132,14 @@ int main(int argc, char *argv[]) {
usage();
/* smoothness */
- else if (argv[fa][1] == 's' || argv[fa][1] == 'S') {
+ else if (argv[fa][1] == 's') {
fa = nfa;
if (na == NULL) usage();
smooth = atof(na);
}
- else if (argv[fa][1] == '2') {
- twopass = 1;
- }
- else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') {
- extra = 1;
+ else if (argv[fa][1] == 'x') {
+ autosm = 1;
}
else
usage();
@@ -225,10 +220,9 @@ int main(int argc, char *argv[]) {
#endif
/* Fit to scattered data */
rss->fit_rspl(rss,
- 0 | (twopass ? RSPL_2PASSSMTH : 0)
- | (extra ? RSPL_EXTRAFIT2 : 0) ,
+ 0 | (autosm ? RSPL_AUTOSMOOTH : 0) ,
test_points, /* Test points */
- pnts, /* Number of test points */
+ pnts, /* Number of test points */
low, high, gres, /* Low, high, resolution of grid */
NULL, NULL, /* Default data scale */
smooth, /* Smoothing */
@@ -245,10 +239,10 @@ int main(int argc, char *argv[]) {
tp.p[0] = x;
rss->interp(rss, &tp);
yy[1][i] = tp.v[0];
- if (yy[1][i] < -0.2)
- yy[1][i] = -0.2;
- else if (yy[1][i] > 1.2)
- yy[1][i] = 1.2;
+ if (yy[1][i] < -20.0)
+ yy[1][i] = -20.0;
+ else if (yy[1][i] > 120.0)
+ yy[1][i] = 120.0;
}
do_plot(xx,yy[0],yy[1],NULL,XRES);
@@ -283,8 +277,7 @@ int main(int argc, char *argv[]) {
gresses[j] = gres[0];
rss->fit_rspl(rss,
- 0 | (twopass ? RSPL_2PASSSMTH : 0)
- | (extra ? RSPL_EXTRAFIT2 : 0) ,
+ 0 | (autosm ? RSPL_AUTOSMOOTH : 0) ,
test_points, /* Test points */
pnts, /* Number of test points */
low, high, gres, /* Low, high, resolution of grid */
@@ -302,10 +295,10 @@ int main(int argc, char *argv[]) {
tp.p[0] = x;
rss->interp(rss, &tp);
yy[1+j][i] = tp.v[0];
- if (yy[1+j][i] < -0.2)
- yy[1+j][i] = -0.2;
- else if (yy[1+j][i] > 1.2)
- yy[1+j][i] = 1.2;
+ if (yy[1+j][i] < -20.0)
+ yy[1+j][i] = -20.0;
+ else if (yy[1+j][i] > 120.0)
+ yy[1+j][i] = 120.0;
}
}
diff --git a/rspl/c1df.c b/rspl/c1df.c
index 5d8dc32..4020d7a 100644
--- a/rspl/c1df.c
+++ b/rspl/c1df.c
@@ -30,8 +30,9 @@
#include "copyright.h"
#include "aconfig.h"
#include "numlib.h"
-#include "plot.h"
#include "rspl.h"
+#include "plot.h"
+#include "ui.h"
double lin();
void usage(void);
diff --git a/rspl/cw1.c b/rspl/cw1.c
index 75a4f55..2a026c8 100644
--- a/rspl/cw1.c
+++ b/rspl/cw1.c
@@ -4,7 +4,7 @@
/************************************************/
/* Discrete regularized spline versions */
-/* Test variable grid spacing in 1D
+/* Test variable grid spacing in 1D */
/* Author: Graeme Gill
* Date: 4/10/95
@@ -28,8 +28,9 @@
#include "copyright.h"
#include "aconfig.h"
#include "numlib.h"
-#include "plot.h"
#include "rspl.h"
+#include "plot.h"
+#include "ui.h"
double lin();
void usage(void);
@@ -230,10 +231,10 @@ int main(int argc, char *argv[]) {
}
}
- printf("Full scale: Black = lin, Red = even, Green = pow2\n");
+ printf("Full scale: Black = ref, 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");
+ printf("Magnified: Black = ref, Red = even, Green = pow2\n");
do_plot6(xx2,yy2[0],yy2[1],yy2[2],NULL,NULL,NULL,XRES);
}
return 0;
diff --git a/rspl/cw3.c b/rspl/cw3.c
index 03c9b81..e26cb7a 100644
--- a/rspl/cw3.c
+++ b/rspl/cw3.c
@@ -4,7 +4,7 @@
/************************************************/
/* Discrete regularized spline versions */
-/* Test variable grid spacing in 3D
+/* Test variable grid spacing in 3D */
/* Author: Graeme Gill
* Date: 28/11/2013
@@ -28,8 +28,9 @@
#include "aconfig.h"
#include "counters.h"
#include "numlib.h"
-#include "plot.h"
#include "rspl.h"
+#include "plot.h"
+#include "ui.h"
double lin();
void usage(void);
diff --git a/rspl/gam.c b/rspl/gam.c
index f8d02d0..e46780f 100644
--- a/rspl/gam.c
+++ b/rspl/gam.c
@@ -4,7 +4,7 @@
* Multi-dimensional regularized spline data structure
*
* Precice gamut surface, gamut pruning, ink limiting and K min/max
- * support routine.s
+ * support routines.
*
* Author: Graeme W. Gill
* Date: 2008/11/21
@@ -45,7 +45,7 @@
#include "sort.h" /* Heap sort */
#include "counters.h" /* Counter macros */
-#include "vrml.h" /* If there is VRML stuff here */
+#include "vrml.h" /* If there is VRML/X3D stuff here */
/* Print a vectors value */
#define DBGVI(text, dim, out, vec, end) \
@@ -879,8 +879,10 @@ for (ss = 0; ss < 2; ss++) { /* Negative side then positive */
rtri *tp;
double col[3], pos[3];
- if ((wrl = new_vrml("gam_diag.wrl", 1)) == NULL)
- error("new_vrml failed\n");
+ if ((wrl = new_vrml("gam_diag", 1, vrml_lab)) == NULL)
+ error("new_vrml failed for 'gam_diag%s\n",vrml_ext());
+
+ wrl->start_line_set(wrl, 0);
/* Display the grid */
EC_INIT(gc);
@@ -906,25 +908,25 @@ for (ss = 0; ss < 2; ss++) { /* Negative side then positive */
pos[2] = (double)gp[2];
if (s->gam.outf != NULL)
s->gam.outf(s->gam.cntxf, pos, pos); /* Apply output lookup */
- wrl->add_vertex(wrl, pos);
+ wrl->add_vertex(wrl, 0, pos);
pos[0] = (double)sp[0];
pos[1] = (double)sp[1];
pos[2] = (double)sp[2];
if (s->gam.outf != NULL)
s->gam.outf(s->gam.cntxf, pos, pos); /* Apply output lookup */
- wrl->add_vertex(wrl, pos);
+ wrl->add_vertex(wrl, 0, pos);
}
DC_INC(cc);
}
EC_INC(gc);
}
- wrl->make_lines(wrl, 2);
- wrl->clear(wrl);
+ wrl->make_lines(wrl, 0, 2);
+ wrl->start_line_set(wrl, 0);
/* Display the current triangles transparently */
if (s->gam.ttop != NULL) {
for (vp = s->gam.vtop; vp != NULL; vp = vp->list) {
- wrl->add_vertex(wrl, vp->v);
+ wrl->add_vertex(wrl, 0, vp->v);
}
/* Set the triangles */
@@ -933,19 +935,19 @@ for (ss = 0; ss < 2; ss++) { /* Negative side then positive */
ix[0] = tp->v[0]->n;
ix[1] = tp->v[1]->n;
ix[2] = tp->v[2]->n;
- wrl->add_triangle(wrl, ix);
+ wrl->add_triangle(wrl, 0, ix);
}
-// wrl->make_triangles(wrl, 0.3, NULL);
- wrl->make_triangles(wrl, 0.0, NULL);
- wrl->clear(wrl);
+// wrl->make_triangles(wrl, 0, 0.3, NULL);
+ wrl->make_triangles(wrl, 0, 0.0, NULL);
+ wrl->start_line_set(wrl, 0);
}
/* Show the active edge as a red cone */
col[0] = 1.0; col[1] = 0.0; col[2] = 0.0; /* Red */
wrl->add_cone(wrl, ep->v[0]->v, ep->v[1]->v, col, 1.0);
-printf("~1 edge v1 = %d = %f %f %f\n", ep->v[0]->gix, ep->v[0]->v[0], ep->v[0]->v[1], ep->v[0]->v[2]);
-printf("~1 edge v2 = %d = %f %f %f\n", ep->v[1]->gix, ep->v[1]->v[0], ep->v[1]->v[1], ep->v[1]->v[2]);
+//printf("~1 edge v1 = %d = %f %f %f\n", ep->v[0]->gix, ep->v[0]->v[0], ep->v[0]->v[1], ep->v[0]->v[2]);
+//printf("~1 edge v2 = %d = %f %f %f\n", ep->v[1]->gix, ep->v[1]->v[0], ep->v[1]->v[1], ep->v[1]->v[2]);
/* Show the candidate verticies */
for (ss = 0; ss < 2; ss++) {
@@ -959,6 +961,7 @@ printf("~1 edge v2 = %d = %f %f %f\n", ep->v[1]->gix, ep->v[1]->v[0], ep->v[1]->
}
}
wrl->del(wrl);
+ printf("Waiting for return after writing gam_diag%s\n",vrml_ext());
getchar();
}
#endif
@@ -1184,7 +1187,7 @@ void rspl_gam_plot(rspl *s, char *name) {
vrml *wrl;
if ((wrl = new_vrml(name, 1, 0)) == NULL)
- error("new_vrml failed\n");
+ error("new_vrml failed for '%s%s'\n",name,vrml_ext());
/* Set the verticies */
for (vp = s->gam.vtop; vp != NULL; vp = vp->list) {
diff --git a/rspl/revbench.c b/rspl/revbench.c
index d0c7566..81e668e 100644
--- a/rspl/revbench.c
+++ b/rspl/revbench.c
@@ -22,6 +22,7 @@
#include "aconfig.h"
#include "rspl.h"
#include "numlib.h"
+//#include "ui.h"
#undef DOLIMIT /* Define to have ink limit */
#define LIMITVAL 2.5 /* Total ink limit sum */
diff --git a/rspl/rspl.c b/rspl/rspl.c
index 23ef6b4..431a7e3 100644
--- a/rspl/rspl.c
+++ b/rspl/rspl.c
@@ -72,6 +72,7 @@ static int re_set_rspl(struct _rspl *s, int flags, void *cbntx,
void (*func)(void *cbntx, double *out, double *in));
static void scan_rspl(struct _rspl *s, int flags, void *cbntx,
void (*func)(void *cbntx, double *out, double *in));
+static int tune_value(struct _rspl *s, co *p);
static void filter_rspl(struct _rspl *s, int flags, void *cbctx,
void (*func)(void *cbntx, float **out, double *in, int cvi));
@@ -172,6 +173,7 @@ printf("!!!! rspl.c using interp_rspl_nl !!!!");
s->set_rspl = set_rspl;
s->scan_rspl = scan_rspl;
s->re_set_rspl = re_set_rspl;
+ s->tune_value = tune_value;
s->opt_rspl = opt_rspl_imp;
s->filter_rspl = filter_rspl;
s->get_in_range = get_in_range;
@@ -1242,6 +1244,153 @@ void (*func)(void *cbntx, double *out, double *in) /* Function to get from */
/* ============================================ */
+/* Tune a single value for simplex interpolation. */
+/* Return 0 on success, 1 if input clipping occured, 2 if output clipping occured */
+static int tune_value (
+struct _rspl *s, /* Pointer to Lut object */
+co *p /* Target value */
+) {
+ int e, di = s->di;
+ int f, fdi = s->fdi;
+ double we[MXDI]; /* Coordinate offset within the grid cell */
+ int si[MXDI]; /* we[] Sort index, [0] = smallest */
+ float *gp; /* Pointer to grid cube base */
+ int rv = 0; /* Register clip */
+
+ /* We are using a simplex (ie. tetrahedral for 3D input) interpolation. */
+
+ DEBLU(("In %s\n", icmPdv(di, p->p)));
+
+ /* Figure out which grid cell the point falls into */
+ {
+ gp = s->g.a; /* Base of grid array */
+ for (e = 0; e < di; e++) {
+ int gres_1 = s->g.res[e]-1;
+ double pe, t;
+ int mi;
+ pe = p->p[e];
+ if (pe < s->g.l[e]) { /* Clip to grid */
+ pe = s->g.l[e];
+ rv = 1;
+ }
+ if (pe > s->g.h[e]) {
+ pe = s->g.h[e];
+ rv = 1;
+ }
+ t = (pe - s->g.l[e])/s->g.w[e];
+ mi = (int)floor(t); /* Grid coordinate */
+ if (mi < 0) /* Limit to valid cube base index range */
+ mi = 0;
+ else if (mi >= gres_1)
+ mi = gres_1-1;
+ gp += mi * s->g.fci[e]; /* Add Index offset for grid cube base in dimen */
+ we[e] = t - (double)mi; /* 1.0 - weight */
+//if (rspldb && di == 3) printf("~1 e = %d, ix = %d, we = %f\n", e, mi, we[e]);
+ }
+ DEBLU(("ix %d, we %s\n", (gp - s->g.a)/s->g.pss, icmPdv(di, p->p)));
+ }
+
+ /* Do selection sort on coordinates */
+ {
+ for (e = 0; e < di; e++)
+ si[e] = e; /* Initial unsorted indexes */
+ for (e = 0; e < (di-1); e++) {
+ double cosn;
+ cosn = we[si[e]]; /* Current smallest value */
+ for (f = e+1; f < di; f++) { /* Check against rest */
+ int tt;
+ tt = si[f];
+ if (cosn > we[tt]) {
+ si[f] = si[e]; /* Exchange */
+ si[e] = tt;
+ cosn = we[tt];
+ }
+ }
+ }
+ }
+ DEBLU(("si[] = %s\n", icmPiv(di, si)));
+
+ /* Now compute the weightings, simplex vertices and output values */
+ {
+ double w, ww = 0.0; /* Current vertex weight, sum of weights squared */
+ double cout[MXDO]; /* Current output value */
+ float *ogp = gp; /* Pointer to grid cube base */
+
+ w = 1.0 - we[si[di-1]]; /* Vertex at base of cell */
+ ww += w * w; /* Sum of weights squared */
+ for (f = 0; f < fdi; f++)
+ cout[f] = w * gp[f];
+
+ DEBLU(("ix %d: w %f * val %s\n", (gp - s->g.a)/s->g.pss, w, icmPfv(fdi,gp)));
+
+ for (e = di-1; e > 0; e--) { /* Middle verticies */
+ w = we[si[e]] - we[si[e-1]];
+ ww += w * w; /* Sum of weights squared */
+ gp += s->g.fci[si[e]]; /* Move to top of cell in next largest dimension */
+ for (f = 0; f < fdi; f++)
+ cout[f] += w * gp[f];
+ DEBLU(("ix %d: w %f * val %s\n", (gp - s->g.a)/s->g.pss, w, icmPfv(fdi,gp)));
+ }
+
+ w = we[si[0]];
+ ww += w * w; /* Sum of weights squared */
+ gp += s->g.fci[si[0]]; /* Far corner from base of cell */
+ for (f = 0; f < fdi; f++)
+ cout[f] += w * gp[f];
+ DEBLU(("ix %d: w %f * val %s\n", (gp - s->g.a)/s->g.pss, w, icmPfv(fdi,gp)));
+ DEBLU(("cout %s\n", icmPdv(fdi, cout)));
+
+ /* We distribute the correction needed in proportion to the */
+ /* interpolation weighting, so the biggest correction is to the */
+ /* closest vertex. */
+ for (f = 0; f < fdi; f++)
+ cout[f] = (p->v[f] - cout[f])/ww; /* Amount to distribute */
+
+ gp = ogp;
+ w = 1.0 - we[si[di-1]]; /* Vertex at base of cell */
+ for (f = 0; f < fdi; f++) {
+ gp[f] += w * cout[f]; /* Apply correction */
+ if (gp[f] < s->g.fmin[f]) {
+ gp[f] = s->g.fmin[f];
+ rv |= 2;
+ } else if (gp[f] > s->g.fmax[f]) {
+ gp[f] = s->g.fmax[f];
+ rv |= 2;
+ }
+ }
+
+ for (e = di-1; e > 0; e--) { /* Middle verticies */
+ w = we[si[e]] - we[si[e-1]];
+ gp += s->g.fci[si[e]]; /* Move to top of cell in next largest dimension */
+ for (f = 0; f < fdi; f++) {
+ gp[f] += w * cout[f]; /* Apply correction */
+ if (gp[f] < s->g.fmin[f]) {
+ gp[f] = s->g.fmin[f];
+ rv |= 2;
+ } else if (gp[f] > s->g.fmax[f]) {
+ gp[f] = s->g.fmax[f];
+ rv |= 2;
+ }
+ }
+ }
+
+ w = we[si[0]];
+ gp += s->g.fci[si[0]]; /* Far corner from base of cell */
+ for (f = 0; f < fdi; f++) {
+ gp[f] += w * cout[f]; /* Apply correction */
+ if (gp[f] < s->g.fmin[f]) {
+ gp[f] = s->g.fmin[f];
+ rv |= 2;
+ } else if (gp[f] > s->g.fmax[f]) {
+ gp[f] = s->g.fmax[f];
+ rv |= 2;
+ }
+ }
+ }
+ return rv;
+}
+
+/* ============================================ */
/* Allow the grid values to be filtered. */
/* For each grid value, provide the input value and */
/* pointers to all the output values in a 3^di grid around */
diff --git a/rspl/rspl.h b/rspl/rspl.h
index 5bccba7..ae0a165 100644
--- a/rspl/rspl.h
+++ b/rspl/rspl.h
@@ -92,7 +92,7 @@ struct _rpnts {
double p[MXRI]; /* Data position [di] */
double v[MXRO]; /* Data value [fdi] */
double k[MXRO]; /* Weight factor (nominally 1.0, less for lower confidence data point) */
- double cv[MXRO]; /* Extra fit corrected v[fdi] */
+// double fe; /* Fit error in output pass (ausm) */
}; typedef struct _rpnts rpnts;
/* Hermite interpolation magic data */
@@ -106,6 +106,14 @@ typedef struct {
#include "rev.h" /* Reverse interpolation defintions */
#include "gam.h" /* Gamut defintions */
+/* Sub-fit itteration information */
+typedef struct {
+ int niters; /* Number of multigrid itterations needed */
+ int **ires; /* Resolution for each itteration and dimension */
+ void **mgtmps[MXRO]; /* Store pointers to re-usable mgtmp when incremental */
+ /* (These don't seem to be used anymore. was incremental removed ?) */
+} it_info;
+
/* Structure for final resolution multi-dimensional regularized spline data */
struct _rspl {
@@ -128,19 +136,19 @@ struct _rspl {
/* Function to set from */
/* Scattered Data point related information */
- int zf; /* Extra fitting flag - Compensate for data fit errors each round */
- int tpsm; /* Two pass smoothing flag (if set to 1). */
- int tpsm2; /* Two pass smoothing 2nd pass flag */
+ int ausm; /* Automatic smoothing enabled flag. */
struct {
int no; /* Number of data points in array */
rpnts *a; /* Array of data points */
- datao vl,vw; /* Data value low/width - used to normalize smoothness values */
+ datao vl, vw; /* Data value low/width - not used */
datao va; /* Data value averages */
+// double fea; /* Fit error average */
} d;
- int niters; /* Number of multigrid itterations needed */
- int **ires; /* Resolution for each itteration and dimension */
- void **mgtmps[MXRO]; /* Store pointers to re-usable mgtmp when incremental */
+
+ it_info ii; /* Main itteration information for final rspl */
+ it_info as_ii; /* Automatic smoothing pre-fit itteration info */
+ it_info asm_ii; /* Automatic smoothing map itteration info */
/* Grid points data */
struct {
@@ -158,7 +166,6 @@ struct _rspl {
/* gres[] entries per dimension. Allows for the possibility of */
/* a non-uniform grid spacing, by adjusting the curvature evaluation */
/* appropriately. */
- double **ccv; /* Curvature compensation array for current outpu (May be NULL) */
int fminmax_valid; /* Min/max/scale cached values valid flag. */
int limitv_cached; /* Flag: Ink limit values have been set in the grid array */
@@ -235,12 +242,8 @@ struct _rspl {
void (*del)(struct _rspl *ss);
/* Combination lags used by various functions */
- /* NOTE that RSPL_2PASSSMTH and RSPL_EXTRAFIT2 are available, but the smoothing */
- /* factors are not setup for them, and they are not sufficiently different from the */
- /* default smoothing to be useful. */
#define RSPL_NOFLAGS 0x0000
-#define RSPL_2PASSSMTH 0x0001 /* Use gaussian filter in 2nd pass to smooth */
-#define RSPL_EXTRAFIT2 0x0002 /* Compensate for data errors each round */
+#define RSPL_AUTOSMOOTH 0x0001 /* Automatically determin local optimal avgdev smoothing */
#define RSPL_SYMDOMAIN 0x0004 /* Maintain symetric smoothness with nonsym. resolution */
#define RSPL_SET_APXLS 0x0020 /* For set_rspl, adjust samples for aproximate least squares */
#define RSPL_FASTREVSETUP 0x0010 /* Do a fast reverse setup at the cost of subsequent speed */
@@ -409,6 +412,13 @@ struct _rspl {
void (*func)(void *cbntx, double *out, double *in) /* Function that gets given values */
);
+ /* Tune a single value. */
+ /* Return 0 on success, 1 if input clipping occured, 2 if output clipping occured */
+ int (*tune_value) (
+ struct _rspl *s, /* Pointer to Lut object */
+ co *p /* Target value */
+ );
+
/* Set values by multi-grid optimisation using the provided function. */
int (*opt_rspl)(
struct _rspl *s,/* this */
diff --git a/rspl/rspl1.h b/rspl/rspl1.h
index d5ea0b9..747abd5 100644
--- a/rspl/rspl1.h
+++ b/rspl/rspl1.h
@@ -44,7 +44,7 @@ typedef struct {
#define RSPL_NOFLAGS 0
struct _rspl {
-
+
/* Private: */
int nig; /* number in interpolation grid */
double gl,gh,gw;/* Interpolation grid scale low, high, grid cell width */
diff --git a/rspl/scat.c b/rspl/scat.c
index ec7e99c..f1535fa 100644
--- a/rspl/scat.c
+++ b/rspl/scat.c
@@ -1,4 +1,6 @@
+// current scat.c with V1.6.3 position curve rspl setup code
+
/*
* Argyll Color Correction System
* Multi-dimensional regularized splines data fitter
@@ -126,12 +128,18 @@
#define DEFAVGDEV 0.5 /* Default average deviation % */
/* 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 ENABLE_2PASSSMTH /* [Define] Enable 2 pass smooth using Gaussian filter */
-#define ENABLE_EXTRAFIT /* [Undef] Enable the extra fit option. Good to combat high smoothness. */
+#undef SMOOTH2 /* [Undef] INCOMPLETE Use 3nd order smoothness rather than curvature. */
+ /* 2nd order is optimal about 2.5 x lower than 3rd order, */
+ /* so an even split between 3rd:2nd would be 1.0 0.4, */
+ /* a 9:1 split would be 0.9 0.04 */
+ /* This also disables the incorrect scaling of smoothness with */
+ /* output range */
+#undef AUTOSM /* [Undef] INCOMPLETE Support auto smoothing using LOOCV */
+
+# define CW2 0.9
+# define CW ((1.0 - CW2) * 0.4)
-#define TWOPASSORDER 2.0 /* Filter order. 2 = Gaussian */
/* Tuning parameters */
#ifdef NEVER
@@ -144,13 +152,12 @@
#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 1.5 /* Multi-grid resolution ratio */
-#undef OVERRLX /* Use over relaxation factor when progress slows (worse accuracy ?) */
+#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. */
#define MAXNI 16 /* Maximum itteration without checking progress */
//#define SMOOTH 0.000100 /* Set nominal smoothing (1.0) */
#define WEAKW 0.1 /* Weak default function nominal effect (1.0) */
-#define ZFCOUNT 1 /* Extra fit repeats */
#else
@@ -166,7 +173,6 @@
#define MAXNI 16 /* [16] Maximum itteration without checking progress */
//#define SMOOTH 0.000100 /* Set nominal smoothing (1.0) */
#define WEAKW 0.1 /* [0.1] Weak default function nominal effect (1.0) */
-#define ZFCOUNT 1 /* [1] Extra fit repeats */
#endif
@@ -200,6 +206,9 @@ struct _mgtmp {
/* Scattered data fit stuff */
struct {
double cw[MXDI]; /* Curvature weight factor */
+#ifdef SMOOTH2
+ double cw2[MXDI]; /* Smoothness weight factor */
+#endif
} sf;
/* Grid points data */
@@ -225,7 +234,6 @@ struct _mgtmp {
/* Equation Solution related (Grid point solutions) */
struct {
- double **ccv; /* [gno][di] Curvature Compensation Values */
double **A; /* A matrix of interpoint weights A[g.no][q.acols] */
int acols; /* A matrix columns needed */
/* Packed indexes run from 0..acols-1 */
@@ -233,17 +241,61 @@ struct _mgtmp {
/* and +/-1 offset in all dimensions, but only the +ve offset */
/* half of the sparse matrix is stored, due to equations */
/* being symetrical. */
- int xcol[HACOMPS+8];/* A array column translation from packed to sparse index */
+#define XCOLPMAX (HACOMPS+8)
+ int xcol[XCOLPMAX]; /* A array column translation from packed to sparse index */
int *ixcol; /* A array column translation from sparse to packed index */
double *b; /* b vector for RHS of simultabeous equation b[g.no] */
double normb; /* normal of b vector */
double *x; /* x solution to A . x = b */
} q;
+#ifdef AUTOSM
+ struct _loocv *as; /* Auto Smooth Setup information */
+#endif
+
}; typedef struct _mgtmp mgtmp;
/* ================================================= */
+
+#ifdef AUTOSM
+
+/* Structure to hold LOOCV equations for multi-grid calculations */
+/* One is created for each resolution. Only used in this file. */
+struct _loocv {
+ mgtmp *m; /* Associated mgtmp */
+
+ int ndcells; /* Number of cells with at least one data point */
+ int *dlist; /* Index of base vertex of cell with at least one data point. */
+ /* ndcells will be filled */
+
+ int *vtx_dlist; /* For base vertex index store the index of the first data point */
+ /* within that cell. -1 for no data */
+ int *dat_dlist; /* Index by data point, store the index of the next data point */
+ /* in the list in the cell. -1 for no more data */
+
+ double *sm; /* smoothness map grid data values in log space, 0.0 for none */
+
+ double **As; /* A matrix of smoothness vertex weights */
+ double *bs; /* b vector for RHS of smoothness equation */
+
+#ifdef SMOOTH2
+ double **As2; /* A matrix of smoothness 2 vertex weights */
+ double *bs2; /* b vector for RHS of smoothness 2 equation */
+#endif
+
+ double **Ad; /* A matrix of data vertex weights */
+ double *bd; /* b vector for RHS of data equation */
+
+ double **AA; /* A matrix of vertex weights */
+ double *bb; /* b vector for RHS of equation */
+ double *xx; /* x vector for solution of equation for a given smoothness */
+
+}; typedef struct _loocv loocv;
+
+#endif /* AUTOSM */
+
+/* ================================================= */
/* Temporary arrays used by cj_line(). We try and avoid */
/* allocating and de-allocating these, and merely expand */
/* them as needed */
@@ -256,15 +308,91 @@ static void init_cj_arrays(cj_arrays *ta);
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 mgtmp *new_mgtmp(rspl *s, int gres[MXDI], double smooth, double avgdev, int f, int issm);
static void free_mgtmp(mgtmp *m);
-static void setup_solve(mgtmp *m);
+static void setup_solve(mgtmp *m, mgtmp *sm);
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);
-static void filter_ccv(rspl *s, double stdev);
-static void init_ccv(mgtmp *m);
-static void comp_extrafit_corr(mgtmp *m);
+static double mgtmp_interp(mgtmp *m, double p[MXDI]);
+#ifdef AUTOSM
+static void setup_sutosmsolve(mgtmp *);
+static void comp_fit_errors(mgtmp *m);
+static void plot_mgtmp1(mgtmp *m);
+#endif /* AUTOSM */
+
+static void set_it_info(
+rspl *s,
+int *tres, /* Single dimension grid resolution for each axis */
+it_info *ii /* it_info to set */
+) {
+ int bres; /* Biggest resolution */
+ int sres; /* Starting resolution */
+ double res;
+ double gratio;
+ int i, e, f;
+
+ bres = 0;
+ for (e = 0; e < s->di; e++) {
+ if (tres[e] > bres)
+ bres = tres[e];
+ }
+
+ /* Figure out how many multigrid steps to use */
+#ifdef SMOOTH2
+ sres = 5; /* Start at minimum grid res of 5 */
+#else
+ sres = 4; /* Start at minimum grid res of 4 */
+#endif
+
+ /* Calculate the resolution scaling ratio and number of itters. */
+ gratio = GRATIO;
+ if (((double)bres/(double)sres) <= gratio) {
+ ii->niters = 2;
+ gratio = (double)bres/(double)sres;
+ } else { /* More than one needed */
+ ii->niters = (int)((log((double)bres) - log((double)sres))/log(gratio) + 0.5);
+ gratio = exp((log((double)bres) - log((double)sres))/(double)ii->niters);
+ ii->niters++;
+ }
+
+ /* Allocate space for resolutions */
+ if ((ii->ires = imatrix(0, ii->niters, 0, s->di)) == NULL)
+ error("rspl: malloc failed - ires[][]");
+
+ /* Fill in the resolution values for each itteration */
+ for (res = (double)sres, i = 0; i < ii->niters; i++) {
+ int ires;
+
+ ires = (int)(res + 0.5);
+ for (e = 0; e < s->di; e++) {
+ if ((ires + 1) >= tres[e]) /* If close enough biger than target res. */
+ ii->ires[i][e] = tres[e];
+ else
+ ii->ires[i][e] = ires;
+ }
+ res *= gratio;
+ }
+
+ /* Assert */
+ for (e = 0; e < s->di; e++) {
+ if (ii->ires[ii->niters-1][e] != tres[e])
+ error("rspl: internal error, final res %d != intended res %d\n",
+ ii->ires[ii->niters-1][e], tres[e]);
+ }
+}
+
+static void freeit_info(
+rspl *s,
+it_info *ii
+) {
+ int i, f;
+
+ if (ii->ires != NULL) {
+ free_imatrix(ii->ires, 0, ii->niters, 0, s->di);
+ ii->ires = NULL;
+ }
+}
+
/* Initialise the regular spline from scattered data */
/* Return non-zero if non-monotonic */
@@ -325,12 +453,7 @@ printf("~1 rspl: flags = 0x%x\n",flags);
if (flags & RSPL_NOVERBOSE) /* Turn off progress messages to stdout */
s->verbose = 0;
-#ifdef ENABLE_2PASSSMTH
- s->tpsm = (flags & RSPL_2PASSSMTH) ? 1 : 0; /* Enable 2 pass smoothing */
-#endif
-#ifdef ENABLE_EXTRAFIT
- s->zf = (flags & RSPL_EXTRAFIT2) ? 2 : 0; /* Enable extra fitting effort */
-#endif
+ s->ausm = (flags & RSPL_AUTOSMOOTH) ? 1 : 0; /* Enable auto smoothing */
s->symdom = (flags & RSPL_SYMDOMAIN) ? 1 : 0; /* Turn on symetric smoothness with gres */
/* Save smoothing factor and Average Deviation */
@@ -485,66 +608,14 @@ printf("~1 rspl: flags = 0x%x\n",flags);
/* Allocate the grid data */
alloc_grid(s);
- /* Zero out the re-usable mgtmps */
- for (f = 0; f < s->fdi; f++) {
- s->mgtmps[f] = NULL;
- }
-
- {
- int sres; /* Starting resolution */
- double res;
- double gratio;
-
- /* Figure out how many multigrid steps to use */
- sres = 4; /* Else start at minimum grid res of 4 */
-
- /* Calculate the resolution scaling ratio and number of itters. */
- gratio = GRATIO;
- if (((double)s->g.bres/(double)sres) <= gratio) {
- s->niters = 2;
- gratio = (double)s->g.bres/(double)sres;
- } else { /* More than one needed */
- s->niters = (int)((log((double)s->g.bres) - log((double)sres))/log(gratio) + 0.5);
- gratio = exp((log((double)s->g.bres) - log((double)sres))/(double)s->niters);
- s->niters++;
- }
-
- /* Allocate space for resolutions and mgtmps pointers */
- if ((s->ires = imatrix(0, s->niters, 0, s->di)) == NULL)
- error("rspl: malloc failed - ires[][]");
-
- for (f = 0; f < s->fdi; f++) {
- if ((s->mgtmps[f] = (void *) calloc(s->niters, sizeof(void *))) == NULL)
- error("rspl: malloc failed - mgtmps[]");
- }
-
- /* Fill in the resolution values for each itteration */
- for (res = (double)sres, i = 0; i < s->niters; i++) {
- int ires;
-
- ires = (int)(res + 0.5);
- for (e = 0; e < s->di; e++) {
- if ((ires + 1) >= s->g.res[e]) /* If close enough biger than target res. */
- s->ires[i][e] = s->g.res[e];
- else
- s->ires[i][e] = ires;
- }
- res *= gratio;
- }
-
- /* Assert */
- for (e = 0; e < s->di; e++) {
- if (s->ires[s->niters-1][e] != s->g.res[e])
- error("rspl: internal error, final res %d != intended res %d\n",
- s->ires[s->niters-1][e], s->g.res[e]);
- }
-
- }
+ /* Setup the number of itterations and resolution for each itteration */
+ set_it_info(s, s->g.res, &s->ii);
/* Do the data point fitting */
return add_rspl_imp(s, 0, d, dtp, dno);
}
+/* Weighting adjustment values */
double adjw[21] = {
7.0896971822529019e-278, 2.7480236142217909e+233, 1.4857837676559724e+166,
1.3997102851752585e-152, 1.3987140593588909e-076, 2.8215833239257504e+243,
@@ -555,6 +626,89 @@ double adjw[21] = {
3.2368131456774088e+262, 6.5639459298208554e+045, 2.0087765219520138e-139
};
+/* Do the fitting for one output plane */
+static mgtmp *
+fit_rspl_plane_imp(
+rspl *s, /* this */
+int f, /* Output plane */
+it_info *ii, /* resolution info */
+double smooth, /* Smoothing factor */
+double avgdev, /* Average deviation to use to set smoothness */
+//mgtmp *sm, /* Optional smoothness map */
+cj_arrays *ta /* Temporary array */
+) {
+ int i, nn; /* Multigreid resolution itteration index */
+ mgtmp *pm = NULL, *m = NULL;
+ mgtmp *psm = NULL, *sm = NULL; /* Smoothing map */
+
+ /* For each resolution (itteration) */
+ for (nn = 0; nn < ii->niters; nn++, pm = m) {
+
+ m = new_mgtmp(s, ii->ires[nn], smooth, avgdev, f, 0);
+
+ // ~~~ want setup solve after creating sm,
+ // but can't setup initial values until after setup_solve
+ // because m->q.x needs to be allocated
+ setup_solve(m, sm);
+
+ if (nn == 0) { /* Make sure we have an initial x[] */
+#ifdef AUTOSM
+ Create sm and set initial values
+~~~~
+ sm = new_mgtmp(s, ii->ires[nn], smooth, avgdev, f, 1);
+ for (i = 0; i < sm->g.no; i++)
+ sm->as.sm[i] = m->sf.cw[e]; /* Value from opt_smooth() */
+ /* ?? how to handle cw2[] etc. ? */
+
+#endif /* AUTOSM */
+
+ for (i = 0; i < m->g.no; i++)
+ m->q.x[i] = s->d.va[f]; /* Start with average data value */
+ } else {
+ init_soln(m, pm); /* Scale from previous resolution */
+ free_mgtmp(pm); /* Free previous grid res solution */
+ pm = NULL;
+
+#ifdef AUTOSM
+ init_soln(sm, psm); /* Scale from previous resolution */
+ free_mgtmp(psm); /* Free previous grid res solution */
+ psm = NULL;
+#endif
+ }
+
+ solve_gres(m, ta,
+#if defined(GRADUATED_TOL)
+ TOL * s->g.res[s->g.brix]/s->ires[nn][s->g.brix],
+#else
+ TOL,
+#endif
+ ii->ires[nn][s->g.brix] >= s->g.res[s->g.brix]); /* Use itterative */
+
+#ifdef DEBUG
+ {
+ int k, gno = m->g.no;
+ double *x = m->q.x; /* x vector for result */
+ printf("Plane %d res %d solution values:\n",f,ii->ires[nn][0]);
+ for (k = 0; k < gno; k++)
+ printf("x[%d] = %f\n",k,x[k]);
+ printf("\n");
+ }
+#endif /* DEBUG */
+
+#ifdef AUTOSM
+ if (nn < (ii->niters-1)) {
+
+ setup_sutosmsolve(sm);
+
+ }
+#endif
+ psm = sm;
+ } /* Next resolution */
+
+ /* return the final resolution mgtmp */
+ return m;
+}
+
/* Do the work of initialising from initial data points. */
/* Return non-zero if non-monotonic */
static int
@@ -578,6 +732,10 @@ add_rspl_imp(
return 0;
}
+#ifdef DEBUG
+ printf("add_rspl_imp: flags 0x%x, dno %d, dtp %d\n",flags,dno,dtp);
+#endif
+
/* Allocate space for points */
/* Allocate the scattered data space */
if ((s->d.a = (rpnts *) malloc(sizeof(rpnts) * dno)) == NULL)
@@ -592,9 +750,9 @@ add_rspl_imp(
for (e = 0; e < s->di; e++)
s->d.a[n].p[e] = dp[i].p[e];
for (f = 0; f < s->fdi; f++) {
- s->d.a[n].cv[f] =
s->d.a[n].v[f] = dp[i].v[f];
s->d.a[n].k[f] = 1.0; /* Assume all data points have same weight */
+// s->d.a[n].fe = 0.0;
}
}
} else if (dtp == 1) { /* Per data point weight */
@@ -605,9 +763,9 @@ add_rspl_imp(
for (e = 0; e < s->di; e++)
s->d.a[n].p[e] = dp[i].p[e];
for (f = 0; f < s->fdi; f++) {
- s->d.a[n].cv[f] =
s->d.a[n].v[f] = dp[i].v[f];
s->d.a[n].k[f] = dp[n].w; /* Weight specified */
+// s->d.a[n].fe = 0.0;
}
}
} else { /* Per data point output weight */
@@ -618,9 +776,9 @@ add_rspl_imp(
for (e = 0; e < s->di; e++)
s->d.a[n].p[e] = dp[i].p[e];
for (f = 0; f < s->fdi; f++) {
- s->d.a[n].cv[f] =
s->d.a[n].v[f] = dp[i].v[f];
s->d.a[n].k[f] = dp[n].w[f]; /* Weight specified */
+// s->d.a[n].fe = 0.0;
}
}
}
@@ -628,89 +786,79 @@ add_rspl_imp(
init_cj_arrays(&ta); /* Zero temporary arrays */
- if (s->verbose && s->zf)
- printf("Doing extra fitting\n");
+ if (s->verbose && s->ausm) {
+#ifdef AUTOSM
+ printf("Doing automatic local smoothing optimization\n");
+#else
+ printf("Automatic local smoothing flag ignored !!!\n");
+#endif
+ }
/* Do fit of grid to data for each output dimension */
for (f = 0; f < fdi; f++) {
- int nn = 0; /* Multigreid resolution itteration index */
- int zfcount = ZFCOUNT; /* Number of extra fit adjustments to do */
- int donezf = 0; /* Count - number of extra fit adjustments done */
float *gp;
mgtmp *m = NULL;
- for (donezf = 0; donezf <= s->zf; donezf++) { /* For each extra fit pass */
+#ifdef NEVER // ~~99 remove this
+ mgtmp *sm = NULL; /* Auto smoothness map */
- for (s->tpsm2 = 0; s->tpsm2 <= s->tpsm; s->tpsm2++) { /* For passes of 2 pass smoothing */
+ /* If auto smoothness, create smoothing map */
+ if (s->ausm) {
+ int res[MXDI];
+ int mxres[5] = { 0, 101, 51, 17, 11 };
+ int smres[5] = { 0, 12, 8, 6, 6 };
- /* For each resolution (itteration) */
- for (nn = 0; nn < s->niters; nn++) {
+ /* Set target resolution for initial fit */
+ for (e = 0; e < s->di; e++) {
+ res[e] = s->g.res[e];
- m = new_mgtmp(s, s->ires[nn], f);
- s->mgtmps[f][nn] = (void *)m;
+ if (res[e] > mxres[s->di])
+ res[e] = mxres[s->di];
+ }
- 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);
+ /* Setup the number of itterations and resolution for each itteration */
+ set_it_info(s, res, &s->as_ii);
- if (nn == 0) { /* Make sure we have an initial x[] */
- for (i = 0; i < m->g.no; i++)
- m->q.x[i] = s->d.va[f]; /* Start with average data value */
- } else {
- init_soln(m, s->mgtmps[f][nn-1]); /* Scale from previous resolution */
+printf("~1 s->smooth = %f, avgdev[f] = %f\n",s->smooth, s->avgdev[f]);
- free_mgtmp(s->mgtmps[f][nn-1]); /* Free previous grid res solution */
- s->mgtmps[f][nn-1] = NULL;
- }
+ /* First pass fit with heavy smoothing */
+ m = fit_rspl_plane_imp(s, f, &s->as_ii, 1.0, 0.1, NULL, &ta);
- solve_gres(m, &ta,
-#if defined(GRADUATED_TOL)
- TOL * s->g.res[s->g.brix]/s->ires[nn][s->g.brix],
-#else
- TOL,
-#endif
- s->ires[nn][s->g.brix] >= s->g.res[s->g.brix]); /* Use itterative */
+printf("Initial high smoothing fit of output %d:\n",f);
+plot_mgtmp1(m);
- } /* Next resolution */
+ /* Compute the fit error values from first pass */
+ comp_fit_errors(m); /* Compute correction to data target values */
- if (s->tpsm && s->tpsm2 == 0) {
- double fstdev; /* Filter standard deviation */
-//printf("~1 setting up second pass smoothing !!!\n");
+ free_mgtmp(m);
- /* Compute the curvature compensation values from */
- /* first pass final resolution */
- comp_ccv(m);
+ /* Set target resolution for smoothness map */
+ for (e = 0; e < s->di; e++)
+ res[e] = smres[s->di];
- if (s->smooth >= 0.0) {
- /* Compute from: no dim, no data points, avgdev & extrafit */
- fstdev = 0.05 * s->smooth;
-fprintf(stderr,"~1 !!! Gaussian smoothing not being computed Using default %f !!!\n",fstdev);
- } else { /* Special used to calibrate table */
- fstdev = -s->smooth;
- }
-//fprintf(stderr,"~1 Gaussian smoothing with fstdev %f !!!\n",fstdev);
- /* Smooth the ccv's */
- filter_ccv(s, fstdev);
- }
- } /* Next two pass smoothing pass */
- if (s->zf)
- comp_extrafit_corr(m); /* Compute correction to data target values */
- } /* Next extra fit pass */
-
- /* Clean up after 2 pass smoothing */
- s->tpsm2 = 0;
- if (s->g.ccv != NULL) {
- free_dmatrix(s->g.ccv, 0, s->g.no-1, 0, s->di-1);
- s->g.ccv = NULL;
+ set_it_info(s, res, &s->asm_ii);
+
+ /* Create smoothness map from fit errors */
+ sm = fit_rspl_plane_imp(s, -1, &s->asm_ii, -50000, 0.0, NULL, &ta);
+// sm = fit_rspl_plane_imp(s, -1, &s->asm_ii, 1000.0, 0.1, NULL, &ta);
+printf("Smoothness map for output %d:\n",f);
+plot_mgtmp1(sm);
}
+#endif /* NEVER */
+
+ /* Fit data for this plane */
+ m = fit_rspl_plane_imp(s, f, &s->ii, s->smooth, s->avgdev[f], /* sm, */ &ta);
+//printf("Final fit for output %d:\n",f);
+//plot_mgtmp1(m);
/* Transfer result in x[] to appropriate grid point value */
for (gp = s->g.a, i = 0; i < s->g.no; gp += s->g.pss, i++)
gp[f] = (float)m->q.x[i];
- free_mgtmp(s->mgtmps[f][nn-1]); /* Free final resolution entry */
- s->mgtmps[f][nn-1] = NULL;
+ free_mgtmp(m); /* Free final resolution entry */
+
+// if (sm != NULL) /* Free smoothing map */
+// free_mgtmp(sm);
} /* Next output channel */
@@ -862,25 +1010,8 @@ init_data(rspl *s) {
/* Free the scattered data allocation */
void
free_data(rspl *s) {
- int i, f;
-
- if (s->ires != NULL) {
- free_imatrix(s->ires, 0, s->niters, 0, s->di);
- s->ires = NULL;
- }
- /* Free up mgtmps */
- for (f = 0; f < s->fdi; f++) {
- if (s->mgtmps[f] != NULL) {
- for (i = 0; i < s->niters; i++) {
- if (s->mgtmps[f][i] != NULL) {
- free_mgtmp(s->mgtmps[f][i]);
- }
- }
- free(s->mgtmps[f]);
- s->mgtmps[f] = NULL;
- }
- }
+ freeit_info(s, &s->ii);
if (s->d.a != NULL) { /* Free up the data point data */
free((void *)s->d.a);
@@ -929,13 +1060,28 @@ free_data(rspl *s) {
/* !!! Possible answer - should be using third order differences */
/* for controlling smoothness, not second order (curvature) ?? */
+/* Should adapt smoothing to noise level in different parts of gamut. */
+/* Fix scaling smoothing to data range bug! */
+
+/* Smoothness tweak */
+static double tweak[21] = {
+ 8.0891733310676571e-263, 1.1269230397087924e+243, 5.5667427967136639e+170,
+ 4.6422059659371074e-072, 4.7573037006103243e-038, 2.2050803446598081e-152,
+ 1.9082109674254010e-094, 1.2362202651281476e+262, 1.8334727652805863e+044,
+ 1.7193993129127580e-139, 8.4028172720870109e-316, 7.7791723264393403e-260,
+ 4.5505694361996285e+198, 1.4450789782663302e+214, 4.8548304485951407e-033,
+ 6.0848773033767158e-153, 2.2014810203887549e+049, 6.0451581453053059e-153,
+ 4.5657997262605343e+233, 1.1415770815909824e+243, 2.0087364177250134e-139
+};
+
+#ifndef SMOOTH2
-// ~~99
static double opt_smooth(
rspl *s,
int di, /* Dimensions */
int ndp, /* Number of data points */
- double ad /* Average sample deviation (proportion of input range) */
+ double ad, /* Average sample deviation (proportion of input range) */
+ int f /* Output chanel */
) {
int i;
double nc; /* Normalised sample count */
@@ -1027,22 +1173,13 @@ static double opt_smooth(
}
};
- /* Smoothness tweak */
- static double tweak[21] = {
- 8.0891733310676571e-263, 1.1269230397087924e+243, 5.5667427967136639e+170,
- 4.6422059659371074e-072, 4.7573037006103243e-038, 2.2050803446598081e-152,
- 1.9082109674254010e-094, 1.2362202651281476e+262, 1.8334727652805863e+044,
- 1.7193993129127580e-139, 8.4028172720870109e-316, 7.7791723264393403e-260,
- 4.5505694361996285e+198, 1.4450789782663302e+214, 4.8548304485951407e-033,
- 6.0848773033767158e-153, 2.2014810203887549e+049, 6.0451581453053059e-153,
- 4.5657997262605343e+233, 1.1415770815909824e+243, 2.0087364177250134e-139
- };
-
/* Real world correction factors go here - */
/* None needed at the moment ? */
double rwf[4] = { 1.0, 1.0, 1.0, 1.0 }; /* Factor for each dimension */
-//printf("~1 opt_smooth called with di = %d, ndp = %d, ad = %f\n",di,ndp,ad);
+#ifdef DEBUG
+ printf("opt_smooth called with di = %d, ndp = %d, ad = %e, f = %d\n",di,ndp,ad,f);
+#endif
if (di < 1)
di = 1;
nc = pow((double)ndp, 1.0/(double)di); /* Normalised sample count */
@@ -1088,13 +1225,17 @@ static double opt_smooth(
}
/* Lookup & interpolate the log smoothness factor */
-//printf("~1 di = %d, ncix = %d, adix = %d\n",di,ncix,adix);
+#ifdef DEBUG
+ printf("di = %d, ncix = %d, adix = %d\n",di,ncix,adix);
+#endif
lsm = smf[di][ncix][adix] * ncw * adw;
lsm += smf[di][ncix][adix+1] * ncw * (1.0 - adw);
lsm += smf[di][ncix+1][adix] * (1.0 - ncw) * adw;
lsm += smf[di][ncix+1][adix+1] * (1.0 - ncw) * (1.0 - adw);
-//printf("~1 lsm = %f\n",lsm);
+#ifdef DEBUG
+ printf("lsm = %f\n",lsm);
+#endif
for (tweakf = 0.0, i = 1; i < 21; i++)
tweakf += tweak[i];
@@ -1105,10 +1246,105 @@ static double opt_smooth(
/* and correct for the real world with a final tweak table */
sm *= rwf[di];
-//printf("Got log smth %f, returning %1.9f from ncix %d, ncw %f, adix %d, adw %f\n", lsm, sm, ncix, ncw, adix, adw);
+#ifdef DEBUG
+ printf("opt_smooth got sm %e before output range adjustment\n",sm);
+#endif
+#ifndef SMOOTH2
+ /* This is incorrect, but is built into the tables of releases. */
+ /* It is one of the things stuffing up XYZ profiles. */
+ /* Remove it after re-tuning with SMOOTH2 */
+ sm *= s->d.vw[f]; /* Scale curvature weight for data range */
+#endif
+#ifdef DEBUG
+ printf("opt_smooth returning sm %e after output range adjustment\n",sm);
+#endif
+
+#ifdef DEBUG
+ printf("Got log smth %f, returning %1.9f from ncix %d, ncw %f, adix %d, adw %f\n", lsm, sm, ncix, ncw, adix, adw);
+#endif
+ return sm;
+}
+
+#else /* Smooth 2 */
+
+/* - - - - - - - - - - - - - - - - - - - - - - - -*/
+/* Smooth2 optimal smoothness calculation. */
+/*
+
+ This is for 3rd order smoothness, and uses a set of fitted
+ equations rather than a table.
+
+*/
+
+static double opt_smooth(
+ rspl *s,
+ int di, /* Dimensions */
+ int ndp, /* Number of data points */
+ double ad, /* Average sample deviation (proportion of input range) */
+ int f /* Output chanel */
+) {
+ int i;
+ double tweakf;
+ struct {
+ double nscale, noffset; /* Number data point scale & offset */
+ double nmax, nmin; /* Number data point max and min values */
+ double dscale, doffset; /* Deviation scale & offset */
+ int l_nid;
+ } params[4] = {
+ { -6.80, 3.50, -2.00, -8.00, 1.50, -0.7 }, /* 1d */
+ { -6.00, -7.50, -1.85, -6.40, 1.60, 0.6 }, /* 2d */
+ { -0.84, -0.36, -1.68, -3.70, 1.75, 1.85 }, /* 3d */
+ { -4.00, 11.20, -0.75, -2.30, 1.55, 2.35 } /* 4d */
+ };
+
+ /* Real world correction factors go here - */
+ /* None needed at the moment ? */
+ double rwf[4] = { 0.0, 0.0, 0.0, 0.0 }; /* Log10 factor to add for each dimension */
+
+ double lndp; /* log10 of ndp */
+ double lad; /* log10 of ad */
+ double lmin; /* log10 of minimum level */
+ double sm, lsm; /* log10 of smoothness needed */
+
+ if (di < 1)
+ di = 1;
+ if (di > 4)
+ di = 4;
+ di--; /* Make di 0..3 */
+ if (ad <= 1e-9)
+ ad = 1e-9;
+
+//printf("~1 opt_smooth2 called with di = %d, nodp = %d, avgdev = %f\n",di,ndp,ad);
+
+ lndp = log10(ndp);
+
+ lmin = lndp * params[di].nscale + params[di].noffset;
+ if (lmin > params[di].nmax)
+ lmin = params[di].nmax;
+ if (lmin < params[di].nmin)
+ lmin = params[di].nmin;
+//printf("lmin = %f from lndp\n",lmin,lndp);
+
+ lad = log10(ad);
+ lsm = log10(pow(10.0, lmin) + pow(10.0, lad * params[di].dscale + params[di].doffset));
+//printf("lsm = %f from lmin %f lad %f dscale %f doff %f\n",lmin,lmin,lad,params[di].dscale,params[di].doffset);
+
+ /* and correct for the real world with a final tweak table */
+ lsm += rwf[di];
+
+ for (tweakf = 0.0, i = 1; i < 21; i++)
+ tweakf += tweak[i];
+ tweakf *= tweak[0];
+
+ sm = pow(10.0, lsm * tweakf);
+
+//printf("Got log smth %f, returning sm %1.9f from di %d, nodp %d, avgdev %f\n", lsm, sm,di+1,ndp,ad);
+
return sm;
}
+#endif /* SMOOTH2 */
+
/* - - - - - - - - - - - - - - - - - - - - - - - -*/
/* Multi-grid temp structure (mgtmp) routines */
@@ -1117,7 +1353,10 @@ static double opt_smooth(
static mgtmp *new_mgtmp(
rspl *s, /* associated rspl */
int gres[MXDI], /* resolution to create */
- int f /* output dimension */
+ double smooth, /* Smoothing factor */
+ double avgdev, /* Average deviation to use to set smoothness */
+ int f, /* output dimension */
+ int issm /* We are creating a smoothness map */
) {
mgtmp *m;
int di = s->di;
@@ -1125,6 +1364,9 @@ static mgtmp *new_mgtmp(
int gno, nigc;
int gres_1[MXDI];
int e, g, n, i;
+#ifdef AUTOSM
+ loocv *as = NULL;
+#endif
/* Allocate a structure */
if ((m = (mgtmp *) calloc(1, sizeof(mgtmp))) == NULL)
@@ -1203,13 +1445,13 @@ static mgtmp *new_mgtmp(
}
/* Compute curvature weighting for matching intermediate resolutions for */
- /* the number of grid points curvature that is accuumulated, as well as the */
+ /* the number of grid points curvature that is accumulated, as well as the */
/* geometric effects of a finer fit to the target surface. */
/* This is all to keep the ratio of sum of smoothness error squared */
/* constant in relationship to the sum of data point error squared. */
for (e = 0; e < di; e++) {
double rsm; /* Resolution smoothness factor */
- double smooth;
+ double smval;
if (s->symdom)
rsm = m->g.res[e]; /* Relative final grid size */
@@ -1221,25 +1463,32 @@ static mgtmp *new_mgtmp(
/* (is ^2 for res. * ^2 with error squared) */
rsm /= nigc; /* Average squared non-smoothness */
- /* 2 pass smoothing */
- if (s->tpsm) {
- double lsm;
+ /* Normal */
+ if (smooth >= 0.0) {
- lsm = -6.0;
- if (s->tpsm2 != 0) /* Two pass smoothing second pass */
- lsm += 2.0; /* Use 100 times the smoothness */
- m->sf.cw[e] = pow(10.0, lsm) * rsm;
+ /* Use auto smoothing map in setup */
+// if (sm != NULL) {
+// smval = 1.0;
+// } else {
- /* Normal */
+ /* Table lookup for optimum smoothing factor */
+ smval = opt_smooth(s, di, s->d.no, avgdev, f);
+//printf("~1 opt_smooth returned %e\n",smval);
+#ifdef SMOOTH2
+ m->sf.cw[e] = CW * smooth * smval * rsm;
+ m->sf.cw2[e] = CW2 * smooth * smval * rsm;
+#else
+ m->sf.cw[e] = smooth * smval * rsm;
+//printf("~1 cw[%d] %f = smooth %f * smval %e * rsm %f\n",e,m->sf.cw[e],smooth,smval,rsm);
+#endif
+ /* Special used to calibrate table */
} else {
-
- if (s->smooth >= 0.0) {
- /* Table lookup for optimum smoothing factor */
- smooth = opt_smooth(s, di, s->d.no, s->avgdev[f]);
- m->sf.cw[e] = s->smooth * smooth * rsm;
- } else { /* Special used to calibrate table */
- m->sf.cw[e] = -s->smooth * rsm;
- }
+#ifdef SMOOTH2
+ m->sf.cw[e] = CW * -smooth * rsm;
+ m->sf.cw2[e] = CW2 * -smooth * rsm;
+#else
+ m->sf.cw[e] = -smooth * rsm;
+#endif
}
}
@@ -1248,6 +1497,33 @@ static mgtmp *new_mgtmp(
/* changes in grid resolution and dimensionality. */
m->wdfw = s->weak * WEAKW/(m->g.no * (double)di);
+#ifdef AUTOSM
+ if (s->ausm) {
+ if ((m->as = (loocv *) calloc(1, sizeof(loocv))) == NULL)
+ error("rspl: malloc failed - loocv");
+
+ as = m->as;
+ as->m = m;
+
+ /* Allocate space for arrays keeping track of */
+ /* cells with data points in them. */
+ if ((as->dlist = ivector(0, dno-1)) == NULL)
+ error("rspl: malloc of vector dlist failed");
+
+ if ((as->vtx_dlist = ivector(0,m->g.no-1)) == NULL)
+ error("rspl: malloc of vector vtx_dlist failed");
+
+ for (i = 0; i < m->g.no; i++)
+ as->vtx_dlist[i] = -1;
+
+ if ((as->dat_dlist = ivector(0, dno-1)) == NULL)
+ error("rspl: malloc of vector dat_dlist failed");
+
+ for (i = 0; i < dno; i++)
+ as->dat_dlist[i] = -1;
+ }
+#endif
+
/* Allocate space for auiliary data point related info */
if ((m->d = (struct mgdat *) calloc(dno, sizeof(struct mgdat))) == NULL)
error("rspl: malloc failed - mgtmp");
@@ -1293,10 +1569,19 @@ static mgtmp *new_mgtmp(
printf("%d: %f\n",e,m->d[n].w[e]);
}
#endif /* DEBUG */
+
+#ifdef AUTOSM
+ if (s->ausm && m->as != NULL) {
+ /* Add data point to per cell list */
+ if (as->vtx_dlist[ix] == -1)
+ as->dlist[as->ndcells++] = ix;
+ as->dat_dlist[n] = as->vtx_dlist[ix];
+ as->vtx_dlist[i] = n;
+ }
+#endif /* AUTOSM */
}
/* Set the solution matricies to unalocated */
- m->q.ccv = NULL;
m->q.A = NULL;
m->q.ixcol = NULL;
m->q.b = NULL;
@@ -1305,6 +1590,21 @@ static mgtmp *new_mgtmp(
return m;
}
+#ifdef AUTOSM
+/* Completely free an loocv */
+static void free_loocv(loocv *as) {
+ mgtmp *m = as->m;
+ rspl *s = m->s;
+
+ free_ivector(as->dlist, 0, s->d.no-1);
+ free_ivector(as->vtx_dlist, 0, m->g.no-1);
+ free_ivector(as->dat_dlist, 0, s->d.no-1);
+
+ free(as);
+}
+
+#endif
+
/* Completely free an mgtmp */
static void free_mgtmp(mgtmp *m) {
int e, di = m->s->di, gno = m->g.no;
@@ -1313,16 +1613,91 @@ static void free_mgtmp(mgtmp *m) {
if (m->g.ipos[e] != NULL)
free(m->g.ipos[e]);
}
- if (m->q.ccv != NULL)
- free_dmatrix(m->q.ccv,0,gno-1,0,di-1);
free_dvector(m->q.x,0,gno-1);
free_dvector(m->q.b,0,gno-1);
free((void *)m->q.ixcol);
free_dmatrix(m->q.A,0,gno-1,0,m->q.acols-1);
free((void *)m->d);
+
+#ifdef AUTOSM
+ if (m->as != NULL)
+ free_loocv(m->as);
+#endif
free((void *)m);
}
+#ifdef NEVER
+
+/* Return the curreb A[][] * x[] value. */
+/* (We use this to determine smoothness factor sensitivity for each */
+/* point with just smoothness weighting factors present) */
+static double smth_err
+(
+ double **A, /* Sparse A[][] matrix */
+ double *x, /* x[] matrix */
+ int gno, /* Total number of unknowns */
+ int acols, /* Use colums in A[][] */
+ int *xcol /* sparse expansion lookup array */
+) {
+ int i, k;
+ double rv;
+
+ /* Compute norm of b - A * x */
+ rv = 0.0;
+ for (i = 0; i < gno; i++) {
+ int k0,k1,k2,k3;
+ double sm = 0.0;
+
+ /* Diagonal and to right in 4's */
+ for (k = 0, k3 = i + xcol[k]; k < acols && k3 < gno; k++, k3 = i + xcol[k]) {
+ sm += A[i][k] * x[k3];
+//printf("i %d: A[%d][%d] %f * x[%d] %f = %f\n", i, i, k, A[i][k], k3, x[k3], A[i][k] * x[k3]);
+ }
+
+ /* Left of diagonal in 4's */
+ /* (We take advantage of the symetry: what would be in the row */
+ /* to the left is repeated in the column above.) */
+ for (k = 1, k3 = i-xcol[k]; k < acols && k3 >= 0; k++, k3 = i-xcol[k]) {
+ sm += A[k3][k] * x[k3];
+//printf("i %d: A[%d][%d] %f * x[%d] %f = %f\n", i, k3, k, A[k3][k], k3, x[k3], A[k3][k] * x[k3]);
+ }
+
+ rv += sm * sm;
+ }
+
+ return rv;
+}
+
+
+/* Print out the smoothness error sensitivity for each data location */
+static void
+print_smsens(mgtmp *m) {
+
+ rspl *s = m->s;
+ int di = s->di;
+ int gno = m->g.no, *gres = m->g.res, *gci = m->g.ci;
+ int i;
+ double **A = m->q.A; /* A matrix of interpoint weights */
+ int acols = m->q.acols; /* A matrix columns needed */
+ int *xcol = m->q.xcol; /* A array column translation from packed to sparse index */
+ double *x = m->q.x; /* x vector for result */
+
+ for (i = 0; i < gno; i++)
+ x[i] = 0.0;
+
+ for (i = 0; i < gno; i++) {
+ double ss;
+
+ x[i] = 1.0;
+ ss = smth_err(A, x, gno, acols, xcol);
+ x[i] = 0.0;
+
+ printf("Smoothness sens %d = %e\n",i,ss);
+ }
+}
+
+#endif /* NEVER */
+
/* Initialise the A[][] and b[] matricies ready to solve, given f */
/* (Can be used to re-initialize an mgtmp for changing curve/extra fit factors) */
/* We are setting up the matrix equation Ax = b to solve, where the aim is */
@@ -1355,7 +1730,8 @@ static void free_mgtmp(mgtmp *m) {
*/
static void setup_solve(
- mgtmp *m /* initialized grid temp structure */
+mgtmp *m, /* initialized grid temp structure */
+mgtmp *sm /* Optional smoothing map for ausm mode */
) {
rspl *s = m->s;
int di = s->di;
@@ -1363,10 +1739,9 @@ static void setup_solve(
int *gres = m->g.res, *gci = m->g.ci;
int f = m->f; /* Output dimensions being worked on */
- double **ccv = m->q.ccv; /* ccv vector for adjusting simultabeous equation */
double **A = m->q.A; /* A matrix of interpoint weights */
int acols = m->q.acols; /* A matrix packed columns needed */
- int *xcol = m->q.xcol; /* A array column translation from packed to sparse index */
+ int *xcol = m->q.xcol; /* A array column translation from packed to sparse index */
int *ixcol = m->q.ixcol; /* A array column translation from sparse to packed index */
double *b = m->q.b; /* b vector for RHS of simultabeous equation */
double *x = m->q.x; /* x vector for LHS of simultabeous equation */
@@ -1374,15 +1749,21 @@ static void setup_solve(
double oawt; /* Overall adjustment weight */
double nbsum; /* normb sum */
-//printf("~1 setup_solve got ccv = 0x%x\n",ccv);
+#ifdef DEBUG
+ printf("setup_solve got sm = %p\n",sm);
+#endif
/* Allocate and init the A array column sparse packing lookup and inverse. */
- /* Note that this only works for a minumum grid resolution of 4. */
+ /* Note that this only works for a minumum grid resolution of 4/5. */
/* The sparse di dimension region allowed for is a +/-1 cube around the point */
/* question, plus +/-2 offsets in axis direction only, */
- /* plus +/-3 offset in axis directions if 2nd order smoothing is defined. */
+ /* plus +/-3 offset in axis directions if SMOOTH2 is defined. */
if (A == NULL) { /* Not been allocated previously */
+#ifdef SMOOTH2
DCOUNT(gc, MXDIDO, di, -3, -3, 4); /* Step through +/- 3 cube offset */
+#else
+ DCOUNT(gc, MXDIDO, di, -2, -2, 3); /* Step through +/- 2 cube offset */
+#endif
int ix; /* Grid point offset in grid points */
acols = 0;
@@ -1405,13 +1786,13 @@ static void setup_solve(
/* or if it has exactly one +/-2 and otherwise 0 */
if ((n3 == 0 && n2 == 0)
|| (n2 == 1 && nz == (di-1))
-#ifdef SMOOTH2
|| (n3 == 1 && nz == (di-1))
-#endif /* SMOOTH2*/
) {
for (ix = 0, k = 0; k < di; k++)
ix += gc[k] * gci[k]; /* Multi-dimension grid offset */
if (ix >= 0) {
+ if (acols >= XCOLPMAX)
+ error("rspl internal: exceeded xcol bounds");
xcol[acols++] = ix; /* We only store half, due to symetry */
}
}
@@ -1459,7 +1840,6 @@ static void setup_solve(
}
/* Stash in the mgtmp */
- m->q.ccv = ccv;
m->q.A = A;
m->q.b = b;
m->q.x = x;
@@ -1475,6 +1855,13 @@ static void setup_solve(
b[i] = 0.0;
}
+ /* Compute adjustment cooeficient */
+ {
+ for (oawt = 0.0, i = 1; i < 21; i++)
+ oawt += wvals[i];
+ oawt *= wvals[0];
+ }
+
/* Production version, without extra edge weight */
/* Accumulate curvature dependent factors to the triangular A matrix. */
@@ -1490,9 +1877,9 @@ static void setup_solve(
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]
+ ki-1 * 1 * 2 * eqn(i-1)
+ ki * -2 * 2 * eqn(i)
+ ki+1 * 1 * 2 * eqn(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,
@@ -1503,26 +1890,14 @@ 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.
*/
+#define V17 /* Enable V1.7 code */
{ /* Setting this up from scratch */
- double dwtw[MXDIDO]; /* Density weight normalizer */
- int wtwc[MXDIDO];
+// double dwtw[MXDIDO]; /* Density weight normalizer = average grid widths */
ECOUNT(gc, MXDIDO, di, 0, gres, 0);
- 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 */
-
+#ifdef NEVER /* We're not using density adjustment */
/* Compute the ipos[] weight normalisation factors */
for (e = 0; e < di; e++) {
if (m->g.ipos[e] == NULL)
@@ -1532,108 +1907,305 @@ static void setup_solve(
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] += w;
+ dwtw[e] += 1.0/w;
}
dwtw[e] /= (gres[e] - 1.0); /* Average weights */
+ dwtw[e] = 1.0/dwtw[e];
//printf("dwtw[%d] = %f\n",e,dwtw[e]);
}
+#endif
EC_INIT(gc);
for (i = 0; i < gno; i++) {
+ double smf = 1.0;
+
+#ifdef AUTOSM
+ /* Lookup smoothing factor map */
+ if (sm != NULL) {
+ double p[MXDI];
+ double avgdev;
+ for (e = 0; e < di; e++)
+ p[e] = gc[e]/(gres[e] - 1.0);
+
+ avgdev = mgtmp_interp(sm, p);
+ smf = opt_smooth(s, di, s->d.no, avgdev, f);
+ }
+#endif /* AUTOSM */
for (e = 0; e < di; e++) { /* For each curvature direction */
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 */
+ double cw = smf * 2.0 * m->sf.cw[e]; /* Overall curvature weight */
+//printf("~1 cw %f = smf %f * 2 * sd.cw[%d] %f\n",cw,smf,e,m->sf.cw[e]);
- /* If at least two above lower edge in this dimension */
/* Add influence on Curvature of cell below */
- if ((gc[e]-2) >= 0) {
+ if ((gc[e]-2) >= 0 && (gc[e]+0) < gres[e]) {
/* double kw = cw * gp[UO_C(e,1)].k; */ /* Cell bellow k value */
double kw = cw;
- dw = w1 = 1.0;
+ 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]);
//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);
+ tt = 0.5 * (w0 + w1); /* Local full normalisation */
+ w0 = tt/w0;
w1 = tt/w1;
-#ifndef EN_CW
- w1 = 1.0;
-#endif
-//printf("[%d][%d] dw = %f, w1 = %f\n",gc[e],i,dw,w1);
+//printf("[%d][%d] w1 = %f\n",gc[e],i,w1);
}
- A[i][ixcol[0]] += dw * w1 * w1 * kw;
- if (ccv != NULL)
- b[i] += kw * (w1) * ccv[i - gci[e]][e]; /* Curvature compensation value */
+
+ A[i][ixcol[0]] += w1 * w1 * kw;
+//printf("A[%d][%d] = %f\n",i,0,A[i][ixcol[0]]);
}
- /* 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 */
double kw = cw;
- dw = w0 = w1 = 1.0;
+ 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);
+ tt = 0.5 * (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);
+//printf("[%d][%d] w0 = %f, w1 = %f\n",gc[e],i,w0,w1);
}
- 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 */
+ A[i][ixcol[0]] += -(w0 + w1) * -(w0 + w1) * kw;
+ A[i][ixcol[gci[e]]] += -(w0 + w1) * w1 * kw * oawt;
+//printf("A[%d][%d] = %f\n",i,0,A[i][ixcol[0]]);
+//printf("A[%d][%d] = %f\n",i,1,A[i][ixcol[gci[e]]]);
}
- /* If at least two below the upper edge in this dimension */
/* Add influence on Curvature of cell above */
- if ((gc[e]+2) < gres[e]) {
+ if ((gc[e]+0) >= 0 && (gc[e]+2) < gres[e]) {
/* double kw = cw * gp[UO_C(e,2)].k; */ /* Cell above k value */
double kw = cw;
- dw = w0 = w1 = 1.0;
+ 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]);
//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);
+ tt = 0.5 * (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);
+//printf("[%d][%d] w0 = %f, w1 = %f\n",gc[e],i,w0,w1);
+ }
+
+ A[i][ixcol[0]] += w0 * w0 * kw;
+ A[i][ixcol[1 * gci[e]]] += w0 * -(w0 + w1) * kw;
+ A[i][ixcol[2 * gci[e]]] += w0 * w1 * kw;
+//printf("A[%d][%d] = %f\n",i,0,A[i][ixcol[0]]);
+//printf("A[%d][%d] = %f\n",i,1,A[i][ixcol[gci[e]]]);
+//printf("A[%d][%d] = %f\n",i,2,A[i][ixcol[2 * gci[e]]]);
+ }
+ }
+ EC_INC(gc);
+ }
+ }
+#ifdef DEBUG
+ printf("After adding 2nd order smoothing:\n");
+ for (i = 0; i < gno; i++) {
+ int *xcol = m->q.xcol;
+ printf("b[%d] = %f\n",i,b[i]);
+ for (k = acols-1; k > 0; k--) {
+ if ((i - xcol[k]) >= 0)
+ printf("A[%d][-%d] = %f\n",i,k,A[i-xcol[k]][k]);
+ }
+ for (k = 0; k < acols && (i + xcol[k]) < gno; k++)
+ printf("A[%d][%d] = %f\n",i,k,A[i][k]);
+ printf("\n");
+ }
+#endif /* DEBUG */
+
+#ifdef SMOOTH2
+ /* Accumulate curvature 2nd order 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] - 3 * u[i] + 3 * u[i+1] - u[i+2] )^2
+ with derivatives wrt each node:
+ ki-1 * 1 * 2 * eqn(i)
+ ki * -3 * 2 * eqn(i)
+ ki+1 * 3 * 2 * eqn(i)
+ ki+2 * 1 * 2 * eqn(i)
+
+ Allowing for scaling of each grid difference by w[i-1], w[i] and w[i+1],
+ 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:
+ where w[i+1] corresponds to the width of cell i+1 to i+2,
+ w' = 1/w
+
+ ki * ( w'[i+1] * (u[i+2] - u[i+1])
+ - 2 * w'[i] * (u[i+1] - u[i])
+ + w'[i-1] * (u[i] - u[i-1]) )^2
+
+ multiply out to group the node values:
+ ki * ( w'[i+1] * u[i+2]
+ - (w'[i+1] + 2 * w'[i]) * u[i+1]
+ + (2 * w'[i] + w'[i-1]) * u[i]
+ + w'[i-1] * u[i-1] )^2
+
+ with derivatives wrt each node:
+ ~~~
+
+ */
+
+ { /* Setting this up from scratch */
+ double dwtw[MXDIDO]; /* Density weight normalizer */
+ ECOUNT(gc, MXDIDO, di, 0, gres, 0);
+ EC_INIT(gc);
+
+ /* 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]);
+ }
+
+ /* We setup the equation to be solved for each grid point. */
+ /* Each grid point participates in foure curvature equations, */
+ /* one centered on the grid line below, one that it's the center of, */
+ /* one centered on the grid line above, and one centered on the */
+ /* grid line two above. The equation setup is for the differential */
+ /* for each of these 2nd order curvature equations to be zero. */
+ for (i = 0; i < gno; i++) {
+ double smf = 1.0;
+
+#ifdef AUTOSM
+ /* Lookup smoothing factor map */
+ if (sm != NULL) {
+ double p[MXDI];
+ double avgdev;
+
+ for (e = 0; e < di; e++)
+ p[e] = gc[e]/(gres[e] - 1.0);
+
+ avgdev = mgtmp_interp(sm, p);
+ smf = opt_smooth(s, di, s->d.no, avgdev, f);
+ }
+#endif /* AUTOSM */
+
+ for (e = 0; e < di; e++) {
+ double w0, w1, w2, tt;
+ double cw = smf * 2.0 * m->sf.cw2[e]; /* Overall curvature weight */
+//printf("gno %d dim %d cw %e\n",i,e,cw);
+
+ /* Add influence on Curvature eqation of cell below */
+ if ((gc[e]-3) >= 0 && (gc[e]+0) < gres[e]) {
+ /* double kw = cw * gp[UO_C(e,1)].k; */ /* Cell bellow k value */
+ double kw = cw;
+ w0 = w1 = w2 = 1.0;
+ if (m->g.ipos[e] != NULL) {
+ w0 = fabs(m->g.ipos[e][gc[e]-2] - m->g.ipos[e][gc[e]-3]);
+ w1 = fabs(m->g.ipos[e][gc[e]-1] - m->g.ipos[e][gc[e]-2]);
+ w2 = fabs(m->g.ipos[e][gc[e]+0] - m->g.ipos[e][gc[e]-1]);
+ tt = 1.0/3.0 * (w0 + w1 + w2);
+ w0 = tt/w0;
+ w1 = tt/w1;
+ w2 = tt/w2;
+ }
+ A[i][ixcol[0]] += w2 * w2 * kw;
+//printf("A[%d][%d] = %f\n",i,0,A[i][ixcol[0]]);
+ }
+ /* Add influence on Curvature of this cell */
+ if ((gc[e]-2) >= 0 && (gc[e]+1) < gres[e]) {
+ /* double kw = cw * gp->k; */ /* This cells k value */
+ double kw = cw;
+ w0 = w1 = w2 = 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]);
+ w2 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]+0]);
+ tt = 1.0/3.0 * (w0 + w1 + w2);
+ w0 = tt/w0;
+ w1 = tt/w1;
+ w2 = tt/w2;
+ }
+ A[i][ixcol[0]] += -(2.0 * w1 + w2) * -(2.0 * w1 + w2) * kw;
+ A[i][ixcol[gci[e]]] += -(2.0 * w1 + w2) * w2 * kw;
+//printf("A[%d][%d] = %f\n",i,0,A[i][ixcol[0]]);
+//printf("A[%d][%d] = %f\n",i,1,A[i][ixcol[gci[e]]]);
+ }
+ /* Add influence on Curvature of cell above */
+ if ((gc[e]-1) >= 0 && (gc[e]+2) < gres[e]) {
+ /* double kw = cw * gp[UO_C(e,2)].k; */ /* Cell above k value */
+ double kw = cw;
+ w0 = w1 = w2 = 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]);
+ w2 = fabs(m->g.ipos[e][gc[e]+2] - m->g.ipos[e][gc[e]+1]);
+ tt = 1.0/3.0 * (w0 + w1 + w2);
+ w0 = tt/w0;
+ w1 = tt/w1;
+ w2 = tt/w2;
}
- 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 */
+ A[i][ixcol[0]] += (w0 + 2.0 * w1) * (w0 + 2.0 * w1) * kw;
+ A[i][ixcol[1 * gci[e]]] += (w0 + 2.0 * w1) * -(2.0 * w1 + w2) * kw;
+ A[i][ixcol[2 * gci[e]]] += (w0 + 2.0 * w1) * w2 * kw * oawt;
+//printf("A[%d][%d] = %f\n",i,0,A[i][ixcol[0]]);
+//printf("A[%d][%d] = %f\n",i,1,A[i][ixcol[gci[e]]]);
+//printf("A[%d][%d] = %f\n",i,2,A[i][ixcol[2 * gci[e]]]);
+ }
+ /* Add influence on Curvature of cell two above */
+ if ((gc[e]+0) >= 0 && (gc[e]+3) < gres[e]) {
+ /* double kw = cw * gp[UO_C(e,3)].k; */ /* Cell two above k value */
+ double kw = cw;
+ w0 = w1 = w2 = 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]);
+ w2 = fabs(m->g.ipos[e][gc[e]+3] - m->g.ipos[e][gc[e]+2]);
+ tt = 1.0/3.0 * (w0 + w1 + w2);
+ w0 = tt/w0;
+ w1 = tt/w1;
+ w2 = tt/w2;
+ }
+ A[i][ixcol[0]] += -w0 * -w0 * kw;
+ A[i][ixcol[1 * gci[e]]] += -w0 * (w0 + 2.0 * w1) * kw;
+ A[i][ixcol[2 * gci[e]]] += -w0 * -(2.0 * w1 + w2) * kw;
+ A[i][ixcol[3 * gci[e]]] += -w0 * w2 * kw;
+//printf("A[%d][%d] = %f\n",i,0,A[i][ixcol[0]]);
+//printf("A[%d][%d] = %f\n",i,1,A[i][ixcol[gci[e]]]);
+//printf("A[%d][%d] = %f\n",i,2,A[i][ixcol[2 * gci[e]]]);
+//printf("A[%d][%d] = %f\n",i,3,A[i][ixcol[3 * gci[e]]]);
}
}
EC_INC(gc);
}
}
+#ifdef DEBUG
+ printf("After adding 3rd order smoothing:\n");
+ for (i = 0; i < gno; i++) {
+ int *xcol = m->q.xcol;
+ printf("b[%d] = %f\n",i,b[i]);
+ for (k = acols-1; k > 0; k--) {
+ if ((i - xcol[k]) >= 0)
+ printf("A[%d][-%d] = %f\n",i,k,A[i-xcol[k]][k]);
+ }
+ for (k = 0; k < acols && (i + xcol[k]) < gno; k++)
+ printf("A[%d][%d] = %f\n",i,k,A[i][k]);
+ printf("\n");
+ }
+#endif /* DEBUG */
#ifdef DEBUG
- printf("After adding curvature equations:\n");
+ printf("After adding 2nd and 3rd order smoothing equations:\n");
for (i = 0; i < gno; i++) {
printf("b[%d] = %f\n",i,b[i]);
for (k = 0; k < acols; k++) {
@@ -1642,6 +2214,7 @@ static void setup_solve(
printf("\n");
}
#endif /* DEBUG */
+#endif /* SMOOTH2 */
nbsum = 0.0; /* Zero sum of b[] squared */
@@ -1649,7 +2222,7 @@ static void setup_solve(
/* weak "data point" exactly at each grid point. */
/* (Note we're not currently doing this in a cache friendly order, */
/* and we're calling the function once for each output component..) */
- if (s->dfunc != NULL) { /* Setting this up from scratch */
+ if (s->dfunc != NULL && f >= 0) { /* Setting this up from scratch */
double iv[MXDI], ov[MXDO];
ECOUNT(gc, MXDIDO, di, 0, gres, 0);
EC_INIT(gc);
@@ -1681,7 +2254,6 @@ static void setup_solve(
printf("\n");
}
#endif /* DEBUG */
-
}
/* Accumulate data point dependent factors */
@@ -1698,8 +2270,8 @@ static void setup_solve(
ai = bp + m->g.hi[j]; /* A matrix index */
w = m->d[n].w[j]; /* Base point grid weight */
- d = 2.0 * s->d.a[n].k[f] * w; /* (2.0, w are derivative factors, k data pnt wgt) */
- tt = d * s->d.a[n].cv[f]; /* Change in (corrected) data component */
+ d = 2.0 * s->d.a[n].k[f] * w; /* (2.0, w are derivtv factors, k data pnt wgt) */
+ tt = d * s->d.a[n].v[f]; /* Change in data component */
nbsum += (2.0 * b[ai] + tt) * tt; /* += (b[ai] + tt)^2 - b[ai]^2 */
b[ai] += tt; /* New data component value */
@@ -1735,300 +2307,18 @@ static void setup_solve(
// exit(0);
}
+#ifdef AUTOSM
-/* Given that we've done a complete fit at the current resolution, */
-/* allocate and compute the curvature error of each grid point and put it in */
-/* s->g.ccv[gno][di] */
-static void comp_ccv(
- mgtmp *m /* Solution to use */
-) {
- rspl *s = m->s;
- int gno = m->g.no, *gres = m->g.res, *gci = m->g.ci;
- int di = s->di;
- double *x = m->q.x; /* Grid solution values */
- int f = m->f; /* Output dimensions being worked on */
- int e, i;
-
- ECOUNT(gc, MXDIDO, di, 0, gres, 0);
- EC_INIT(gc);
-
- if (s->g.ccv == NULL) {
- if ((s->g.ccv = dmatrixz(0, gno-1, 0, di-1)) == NULL) {
- error("Malloc of ccv[] failed with [%d][%d]",di,gno);
- }
- }
-
- for (i = 0; i < gno; i++) {
- for (e = 0; e < di; e++) { /* For each curvature direction */
- double w0, w1, tt;
-
- s->g.ccv[i][e] = 0.0; /* Default value */
-
- /* If not one from upper or lower edge in this dimension */
- if ((gc[e]-1) >= 0 && (gc[e]+1) < gres[e]) {
- /* double kw = cw * gp->k; */ /* This cells k value */
- 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;
- }
- s->g.ccv[i][e] += w0 * x[i - gci[e]];
- s->g.ccv[i][e] += -(w0 + w1) * x[i];
- s->g.ccv[i][e] += w1 * x[i + gci[e]];
- }
-//printf("~1 computing ccv for node %d is %f\n",i,s->g.ccv[i][0]);
- }
- EC_INC(gc);
- }
-}
-
-/* Down sample the curvature compensation values in s->g.ccv to */
-/* a given solution. Allocate the m->q.ccv if necessary. */
-static void init_ccv(
- mgtmp *m /* Destination */
-) {
- rspl *s = m->s;
- int f = m->f; /* Output dimensions being worked on */
- int di = s->di;
- int gno = m->g.no;
- int gres1_1[MXDI]; /* Destination */
- int gres2_1[MXDI]; /* Source */
- double scale[MXDI]; /* ccv scale factor */
- int e, n;
- ECOUNT(gc, MXDIDO, di, 0, m->g.res, 0); /* Counter for output points */
-
- for (e = 0; e < di; e++) {
- gres1_1[e] = m->g.res[e]-1;
- gres2_1[e] = s->g.res[e]-1;
- }
-
- if (m->q.ccv == NULL) {
- if ((m->q.ccv = dmatrixz(0, gno-1, 0, di-1)) == NULL) {
- error("Malloc of ccv[] failed with [%d][%d]",di,gno);
- }
- }
-
- /* Compute the scale factor to compensate for the grid resolution */
- /* effect on the grid difference values. */
- for (e = 0; e < di; e++) {
- double rsm_s, rsm_d;
-
- if (s->symdom) { /* Relative final grid size */
- rsm_s = s->g.res[e];
- rsm_d = m->g.res[e];
- } else { /* Relative mean final grid size */
- rsm_s = s->g.mres;
- rsm_d = m->g.mres;
- }
-
- rsm_s = pow((rsm_s-1.0), 2.0); /* Geometric resolution factor for smooth surfaces */
- rsm_d = pow((rsm_d-1.0), 2.0); /* (It's ^2 rather than ^4 as it's is before squaring) */
-
- scale[e] = rsm_s/rsm_d;
- }
-
- /* Point sampling is probably not the ideal way of down sampling, */
- /* but it's easy, and won't be too bad if the s->g.ccv has been */
- /* low pass filtered. */
-
- /* For all grid ccv's */
- EC_INIT(gc);
- for (n = 0; n < gno; n++) {
- double we[MXRI]; /* 1.0 - Weight in each dimension */
- double gw[POW2MXRI]; /* weight for each grid cube corner */
- int ix; /* Index of source ccv grid cube base */
-
- /* Figure out which grid cell the point falls into */
- {
- double t;
- int mi;
- ix = 0;
- for (e = 0; e < di; e++) {
- t = ((double)gc[e]/(double)gres1_1[e]) * (double)gres2_1[e];
- mi = (int)floor(t); /* Grid coordinate */
- if (mi < 0) /* Limit to valid cube base index range */
- mi = 0;
- else if (mi >= gres2_1[e])
- mi = gres2_1[e]-1;
- ix += mi * s->g.ci[e]; /* Add Index offset for grid cube base in dimen */
- we[e] = t - (double)mi; /* 1.0 - weight */
- }
- }
-
- /* Compute corner weights needed for interpolation */
- {
- int i, g;
- gw[0] = 1.0;
- for (e = 0, g = 1; e < di; g *= 2, e++) {
- for (i = 0; i < g; i++) {
- gw[g+i] = gw[i] * we[e];
- gw[i] *= (1.0 - we[e]);
- }
- }
- }
-
- /* Compute the output values */
- {
- int i;
- for (e = 0; e < di; e++)
- m->q.ccv[n][e] = 0.0; /* Zero output value */
- for (i = 0; i < (1 << di); i++) { /* For all corners of cube */
- int oix = ix + s->g.hi[i];
- for (e = 0; e < di; e++)
- m->q.ccv[n][e] += gw[i] * s->g.ccv[oix][e];
- }
- /* Rescale curvature for grid spacing */
- for (e = 0; e < di; e++)
- m->q.ccv[n][e] *= scale[e];
-//printf("~1 downsampling ccv for node %d is %f\n",n,m->q.ccv[n][0]);
- }
- EC_INC(gc);
- }
-}
-
-/* Apply a gaussian filter to the curvature compensation values */
-/* s->g.ccv[gno][di], to apply smoothing. */
-static void filter_ccv(
- rspl *s,
- double stdev /* Standard deviation diameter of filter (in input) */
- /* 1.0 = grid width */
-) {
- int k, e, ee, i, j, di = s->di;
- int gno = s->g.no, *gres = s->g.res, *gci = s->g.ci;
- double *_fkern[MXDI], *fkern[MXDI]; /* Filter kernels */
- int kmin[MXDI], kmax[MXDI]; /* Kernel index range (inclusive) */
- double *_row, *row; /* Extended copy of each row processed */
-
-//printf("Doing filter stdev %f\n",stdev);
-
-//printf("~1 bres = %d, index %d to %d\n",s->g.bres,-s->g.bres+1,s->g.bres+s->g.bres-1);
- if ((_row = (double *) malloc(sizeof(double) * (s->g.bres * 3 - 2))) == NULL)
- error("rspl malloc failed - ccv row copy");
- row = _row + s->g.bres-1; /* Allow +/- gres-1 */
-
- /* Compute the kernel weightings for the given stdev */
- for (ee = 0; ee < di; ee++) { /* For each dimension direction */
- int cres; /* Current res */
- double k1, k2, tot;
-
-//printf("Filter along dim %d:\n",ee);
- if ((_fkern[ee] = (double *) malloc(sizeof(double) * (gres[ee] * 2 - 1))) == NULL)
- error("rspl malloc failed - ccv filter kernel");
- fkern[ee] = _fkern[ee] + gres[ee]-1; /* node of interest at center */
-
- /* Take gaussian constants out of the loop */
- k2 = 1.0 / (2.0 * pow(fabs(stdev), TWOPASSORDER));
- k1 = k2 / 3.1415926;
-
- /* Comute the range needed */
- if (s->symdom) {
- cres = gres[ee];
- } else {
- cres = s->g.mres;
- }
- kmin[ee] = (int)floor(-5.0 * stdev * (cres-1.0));
- kmax[ee] = (int)ceil(5.0 * stdev * (cres-1.0));
-
- if (kmin[ee] < (-gres[ee]+1))
- kmin[ee] = -gres[ee]+1;
- else if (kmin[ee] > -1)
- kmin[ee] = -1;
- if (kmax[ee] > (gres[ee]-1))
- kmax[ee] = gres[ee]-1;
- else if (kmax[ee] < 1)
- kmax[ee] = 1;
-//printf("kmin = %d, kmax = %d\n",kmin[ee], kmax[ee]);
-
- for (tot = 0.0, i = kmin[ee]; i <= kmax[ee]; i++) {
- double fi = (double)i;
-
- /* Do a discrete integration of the gassian function */
- /* to compute discrete weightings */
- fkern[ee][i] = 0.0;
- for (k = -4; k < 5; k++) {
- double oset = (fi + k/9.0)/(cres-1.0);
- double val;
-
- val = k1 * exp(-k2 * pow(fabs(oset), TWOPASSORDER));
- fkern[ee][i] += val;
- tot += val;
- }
- }
- /* Normalize the sum */
- for (tot = 1.0/tot, i = kmin[ee]; i <= kmax[ee]; i++)
- fkern[ee][i] *= tot;
-//printf("Filter cooefs:\n");
-//for (i = kmin[ee]; i <= kmax[ee]; i++)
-//printf("%d: %e\n",i,fkern[ee][i]);
- }
-
- for (k = 0; k < di; k++) { /* For each curvature direction */
- for (ee = 0; ee < di; ee++) { /* For each dimension direction */
- int tgres[MXDI-1];
-
-//printf("~1 Filtering curv dir %d, dim dir %d\n",k,ee);
- /* Setup counters for scanning through all other dimensions */
- for (j = e = 0; e < di; e++) {
- if (e == ee)
- continue;
- tgres[j++] = gres[e];
- }
- /* For each row of this dimension */
- {
- ECOUNT(gc, MXDIDO-1, di-1, 0, tgres, 0); /* Count other dimensions */
-
- EC_INIT(gc);
- for (; di <= 1 || !EC_DONE(gc);) {
- int ix;
-
- /* Compute index of start of row */
- for (ix = j = e = 0; e < di; e++) {
- if (e == ee)
- continue;
- ix += gc[j++] * gci[e];
- }
+~~~~9999
- /* Copy row to temporary array, and expand */
- /* edge values by mirroring them. */
- for (i = 0; i < gres[ee]; i++)
- row[i] = s->g.ccv[ix + i * gci[ee]][k];
- for (i = kmin[ee]; i < 0; i++)
- row[i] = 2.0 * row[0] - row[-i]; /* Mirror the value */
- for (i = gres[ee]-1 + kmax[ee]; i > (gres[ee]-1); i--)
- row[i] = 2.0 * row[gres[ee]-1] - row[gres[ee]-1-i]; /* Mirror the value */
-//printf("~1 Row = \n");
-//for (i = kmin[ee]; i <= (gres[ee]-1 + kmax[ee]); i++)
-//printf("%d: %f\n",i,row[i]);
-
- /* Apply the 1D convolution to the temporary array */
- /* to produce the filtered values. */
- for (i = 0; i < gres[ee]; i++) {
- double fv;
-
- for (fv = 0.0, j = kmin[ee]; j <= kmax[ee]; j++)
- fv += fkern[ee][j] * row[i + j];
- s->g.ccv[ix + i * gci[ee]][k] = fv;
- }
- if (di <= 1)
- break;
- EC_INC(gc);
- }
- }
- }
- }
+#endif
- for (ee = 0; ee < di; ee++)
- free(_fkern[ee] );
- free(_row);
-}
+#ifdef AUTOSM
/* Given that we've done a complete fit at the current resolution, */
-/* compute the error of each data point, and then compute */
-/* a correction factor .cv[] for each point from this. */
-static void comp_extrafit_corr(
+/* compute the error of each data point. */
+/* This is used to compute a smoothness factor map */
+static void comp_fit_errors(
mgtmp *m /* Current resolution mgtmp */
) {
rspl *s = m->s;
@@ -2037,6 +2327,7 @@ static void comp_extrafit_corr(
int di = s->di;
double *x = m->q.x; /* Grid solution values */
int f = m->f; /* Output dimensions being worked on */
+ double fea = 0.0; /* Average value */
/* Compute error for each data point */
for (n = 0; n < dno; n++) {
@@ -2047,34 +2338,74 @@ static void comp_extrafit_corr(
double gain = 1.0;
/* Compute the interpolated grid value for this data point */
- for (val = 0.0, j = 0; j < (1 << di); j++) { /* Binary sequence */
+ for (val = 0.0, j = 0; j < (1 << di); j++) /* Binary sequence */
val += m->d[n].w[j] * x[bp + m->g.hi[j]];
- }
err = s->d.a[n].v[f] - val;
+ err *= 0.8;
+// s->d.a[n].fe = fabs(err);
+//printf("~1 data %d fe = %f\n",n,s->d.a[n].fe);
+ fea += s->d.a[n].fe;
+ }
+ fea /= (double)dno;
-#ifdef NEVER
- /* Compute gain from previous move */
- if (fabs(s->d.a[n].pe[f]) > 0.001) {
- gain = (val - s->d.a[n].pv[f])/s->d.a[n].pe[f];
- if (gain < 0.2)
- gain = 0.2;
- else if (gain > 5.0)
- gain = 5.0;
- gain = pow(gain, 0.6);
- } else {
- gain = 1.0;
+// s->d.fea = fea; /* Average fit error */
+}
+
+#endif /* AUTOSM */
+
+/* Return an interpolayed value */
+static double mgtmp_interp(
+ mgtmp *m,
+ double p[MXDI] /* Input coord in normalised grid forms */
+) {
+ rspl *s = m->s;
+ int di = s->di;
+ int e, n;
+ double we[MXRI]; /* 1.0 - Weight in each dimension */
+ double gw[POW2MXRI]; /* weight for each grid cube corner */
+ double *gp; /* Pointer to x2[] grid cube base */
+ double val;
+
+ /* Figure out which grid cell the point falls into */
+ {
+ double t;
+ int mi;
+ gp = m->q.x; /* Base of solution array */
+ for (e = 0; e < di; e++) {
+ t = (double)p[e] * (m->g.res[e] - 1.0);
+ mi = (int)floor(t); /* Grid coordinate */
+ if (mi < 0) /* Limit to valid cube base index range */
+ mi = 0;
+ else if (mi >= (m->g.res[e] - 1))
+ mi = m->g.res[e] - 2;
+ gp += mi * m->g.ci[e]; /* Add Index offset for grid cube base in dimen */
+ we[e] = t - (double)mi; /* 1.0 - weight */
}
-#endif
- /* Correct the target data point value by the error */
- s->d.a[n].cv[f] += err / gain;
-
-//printf("~1 Data point %d, v = %f, cv = %f, change = %f\n",n,s->d.a[n].v[f],s->d.a[n].cv[f],-val);
-//printf("~1 Data point %d, pe = %f, change = %f, gain = %f\n",n,s->d.a[n].pe[f],val - s->d.a[n].pv[f],gain);
-//printf("~1 Data point %d err = %f, target %f, was %f, now %f\n",n,err,s->d.a[n].v[f],val,s->d.a[n].cv[f]);
-// s->d.a[n].pe[f] = err / gain;
-// s->d.a[n].pv[f] = val;
}
+
+ /* Compute corner weights needed for interpolation */
+ {
+ int i, g;
+ gw[0] = 1.0;
+ for (e = 0, g = 1; e < di; g *= 2, e++) {
+ for (i = 0; i < g; i++) {
+ gw[g+i] = gw[i] * we[e];
+ gw[i] *= (1.0 - we[e]);
+ }
+ }
+ }
+
+ /* Compute the output values */
+ {
+ int i;
+ val = 0.0; /* Zero output value */
+ for (i = 0; i < (1 << di); i++) { /* For all corners of cube */
+ val += gw[i] * gp[m->g.hi[i]];
+ }
+ }
+
+ return val;
}
/* Transfer a solution from one mgtmp to another */
@@ -2086,64 +2417,48 @@ static void init_soln(
rspl *s = m1->s;
int di = s->di;
int gno = m1->g.no;
- int gres1_1[MXDI];
- int gres2_1[MXDI];
int e, n;
ECOUNT(gc, MXDIDO, di, 0, m1->g.res, 0); /* Counter for output points */
- for (e = 0; e < di; e++) {
- gres1_1[e] = m1->g.res[e]-1;
- gres2_1[e] = m2->g.res[e]-1;
- }
-
/* For all output grid points */
EC_INIT(gc);
for (n = 0; n < gno; n++) {
- double we[MXRI]; /* 1.0 - Weight in each dimension */
- double gw[POW2MXRI]; /* weight for each grid cube corner */
- double *gp; /* Pointer to x2[] grid cube base */
-
- /* Figure out which grid cell the point falls into */
- {
- double t;
- int mi;
- gp = m2->q.x; /* Base of solution array */
- for (e = 0; e < di; e++) {
- t = ((double)gc[e]/(double)gres1_1[e]) * (double)gres2_1[e];
- mi = (int)floor(t); /* Grid coordinate */
- if (mi < 0) /* Limit to valid cube base index range */
- mi = 0;
- else if (mi >= gres2_1[e])
- mi = gres2_1[e]-1;
- gp += mi * m2->g.ci[e]; /* Add Index offset for grid cube base in dimen */
- we[e] = t - (double)mi; /* 1.0 - weight */
- }
- }
+ double p[MXRI]; /* Grid relative location */
- /* Compute corner weights needed for interpolation */
- {
- int i, g;
- gw[0] = 1.0;
- for (e = 0, g = 1; e < di; g *= 2, e++) {
- for (i = 0; i < g; i++) {
- gw[g+i] = gw[i] * we[e];
- gw[i] *= (1.0 - we[e]);
- }
- }
- }
+ for (e = 0; e < di; e++)
+ p[e] = (double)gc[e]/(m1->g.res[e] - 1.0);
+
+ m1->q.x[n] = mgtmp_interp(m2, p);
- /* Compute the output values */
- {
- int i;
- m1->q.x[n] = 0.0; /* Zero output value */
- for (i = 0; i < (1 << di); i++) { /* For all corners of cube */
- m1->q.x[n] += gw[i] * gp[m2->g.hi[i]];
- }
- }
EC_INC(gc);
}
}
+
+#ifdef AUTOSM
+
+#ifndef NEVER // Debug
+
+/* Plot the 0'th dimension response */
+void plot_mgtmp1(mgtmp *m) {
+ int i;
+ double xx[100];
+ double yy[100];
+ double p[MXDI];
+
+ for (i = 0; i < m->s->di; i++)
+ p[i] = 0.0;
+
+ for (i = 0; i < 100; i++) {
+ xx[i] = p[0] = (double)i/99.0;
+ yy[i] = mgtmp_interp(m, p);
+ }
+ do_plot(xx, yy, NULL, NULL, 100);
+}
+
+#endif /* NEVER */
+#endif /* AUTOSM */
+
/* - - - - - - - - - - - - - - - - - - - -*/
static double one_itter1(cj_arrays *ta, double **A, double *x, double *b, double normb,
@@ -2234,7 +2549,7 @@ solve_gres(mgtmp *m, cj_arrays *ta, double tol, int final)
/* Note that we process the A[][] sparse columns in compact form */
#ifdef DEBUG_PROGRESS
- printf("Target tol = %f\n",tol);
+ printf("Target tol = %e\n",tol);
#endif
/* If the number of point is small, or it is just one */
/* dimensional, solve it more directly. */
@@ -2417,7 +2732,7 @@ one_itter2(
sm += A[i][k] * x[k3];
/* Left of diagonal in 4's */
- /* (We take advantage of the symetry: what would be in the row */
+ /* (We take advantage of A[][] symetry: what would be in the row */
/* to the left is repeated in the column above.) */
for (k = 1, k3 = i-xcol[k+3]; (k+3) < acols && k3 >= 0; k += 4, k3 = i-xcol[k+3]) {
k0 = i-xcol[k+0];
@@ -2432,6 +2747,7 @@ one_itter2(
for (k3 = i-xcol[k]; k < acols && k3 >= 0; k++, k3 = i-xcol[k])
sm += A[k3][k] * x[k3];
+ /* Compute x value that solves equation just for this point */
// x[i] = (b[i] - sm)/A[i][0];
x[i] += ovsh * ((b[i] - sm)/A[i][0] - x[i]);
@@ -2627,8 +2943,10 @@ cj_line(
sm += A[ii][k] * x[k3];
/* Left of diagonal in 4's */
- /* (We take advantage of the symetry: what would be in the row */
- /* to the left is repeated in the column above.) */
+ /* (We take advantage of the symetry around the diagonal: what would be in the row */
+ /* to the left is repeated in the column above, so for an unsparse matrix */
+ /* we simply swapt the row and column index, for sparse we use the mirror column */
+ /* offset (ie. to the right side) and subtract the column offset from the row.) */
for (k = 1, k3 = ii-xcol[k+3]; (k+3) < acols && k3 >= 0; k += 4, k3 = ii-xcol[k+3]) {
k0 = ii-xcol[k+0];
k1 = ii-xcol[k+1];
@@ -2682,6 +3000,7 @@ cj_line(
max_it = 0;
}
+ /* Improve the solution */
for (it = 1; it <= max_it; it++) {
/* Aproximately solve for z[] given r[], */
diff --git a/rspl/smtmpp.c b/rspl/smtmpp.c
index ed7bac7..1264d4b 100644
--- a/rspl/smtmpp.c
+++ b/rspl/smtmpp.c
@@ -28,9 +28,10 @@
#include "rspl.h"
#include "numlib.h"
#include "xicc.h"
-#include "plot.h"
#include "rspl_imp.h"
#include "counters.h" /* Counter macros */
+#include "plot.h"
+#include "ui.h"
#ifdef DEBUG
#define DBG(xxxx) printf xxxx
@@ -67,6 +68,7 @@ refconv *rco,
double *out,
double *in) {
rco->mppo->lookup(rco->mppo, out, in);
+//printf("~1 3/4D %f %f %f -> %f %f %f\n", in[0], in[1], in[2], out[0], out[1], out[2]);
}
/* Do a 1d emulation */
@@ -129,6 +131,8 @@ static int set_refconv(
s->di = di;
s->ix = ix;
+ s->lookup = NULL;
+
if (di == s->pdi) {
if (ix == 0) {
s->lookup = refconv_default;
@@ -185,9 +189,7 @@ static void do_test(
int ntps, /* Number of sample points */
double noise, /* Sample point noise volume */
int unif, /* nz for uniform noise, else normal */
- double smooth, /* Smoothness to test */
- int twopass, /* Two pass flag */
- int extra /* Extra fit flag */
+ double smooth /* Smoothness to test */
);
/* Compute smoothness of function */
@@ -213,7 +215,7 @@ static double best(int n, double *rmse, double *smv) {
int i, bi;
rspl *curve;
co *tps = NULL;
- int ns = 500; /* Number of steps to search */
+ int ns = 2000; /* Number of samples */
datai low,high;
int gres[1];
datai dlow,dhigh;
@@ -248,7 +250,7 @@ static double best(int n, double *rmse, double *smv) {
n, /* Number of test points */
NULL, NULL, gres, /* Low, high, resolution of grid */
NULL, NULL, /* Default data scale */
- -0.000005, /* Underlying smoothing */
+ -0.0007, /* Underlying smoothing */
avgdev, /* Average deviation */
NULL); /* iwidth */
@@ -262,30 +264,13 @@ static double best(int n, double *rmse, double *smv) {
printf("Point %d at %f, should be %f is %f\n",i,log10(smv[i]),rmse[i],tp.v[0]);
}
-
-#define TPRES 100
- /* Plot the result */
- {
- double xx[TPRES], yy[TPRES];
-
- for (i = 0; i < TPRES; i++) {
- co tp;
- double vi = i/(TPRES-1.0);
-
- tp.p[0] = log10(smv[0]) + (log10(smv[n-1]) - log10(smv[0])) * vi;
- curve->interp(curve, &tp);
- xx[i] = tp.p[0];
- yy[i] = tp.v[0];
- }
- do_plot(xx,yy,NULL,NULL,TPRES);
- }
#endif
/* Choose a solution */
-
- /* First find the very lowest error */
brmse = 1e38;
- for (i = 0; i < ns ; i++) {
+
+ /* Find lowest rms error point */
+ for (i = ns-1; i >= 0; i--) {
co tp;
double vi;
@@ -300,15 +285,8 @@ static double best(int n, double *rmse, double *smv) {
}
}
-//printf("~1 located minimum at %f err %f\n",pow(10.0, blsmv), brmse);
-
- /* Then locate the larger smoothness value that */
- /* gives a slightly higher error. */
- if ((brmse * 1.1) < (brmse + 1.0))
- brmse *= 1.1; /* + 10% */
- else
- brmse += 1.0; /* or 1 delta E */
- for (i = bi; i < ns ; i++) {
+ /* Then increase smoothness until fit error is 1% higher */
+ for (i = bi+1; i < ns; i++) {
co tp;
double vi;
@@ -316,23 +294,41 @@ static double best(int n, double *rmse, double *smv) {
tp.p[0] = log10(smv[0]) + (log10(smv[n-1]) - log10(smv[0])) * vi;
curve->interp(curve, &tp);
- if (tp.v[0] >= brmse) {
+ if (tp.v[0] >= (1.01 * brmse)) {
blsmv = tp.p[0];
brmse = tp.v[0];
break;
}
}
-//printf("~1 located minimum + 20%% at %f err %f\n",pow(10.0, blsmv), brmse);
-
rv = pow(10.0, blsmv);
+
+#ifdef NEVER
+
+#define TPRES 100
+ /* Plot the result */
+ {
+ double xx[TPRES], yy[TPRES];
+
+ for (i = 0; i < TPRES; i++) {
+ co tp;
+ double vi = i/(TPRES-1.0);
+
+ tp.p[0] = log10(smv[0]) + (log10(smv[n-1]) - log10(smv[0])) * vi;
+ curve->interp(curve, &tp);
+ xx[i] = tp.p[0];
+ yy[i] = tp.v[0];
+ }
+ printf("Best at %f\n",blsmv);
+ do_plot(xx,yy,NULL,NULL,TPRES);
+ }
+#endif
+
return rv;
}
/* ---------------------------------------------------------------------- */
/* Explore ideal smoothness change with test point number and noise volume */
-static void do_series_1(refconv* rco, int tdi, int unif, int tntps, int tnlev, int twopass, int extra) {
- int verb = 0;
- int plot = 0;
+static void do_series_1(int verb, int plot, refconv* rco, int tdi, int unif, int tntps, int tnlev) {
int sdi = 1, edi = 4, di;
int res = 0;
int ntps = 0;
@@ -350,111 +346,118 @@ static void do_series_1(refconv* rco, int tdi, int unif, int tntps, int tnlev, i
};
/* Set of smoothnesses to explore */
- double smset[20] = {
- -0.0000001,
- -0.0000010,
- -0.0000050,
- -0.0000100,
- -0.0000500,
- -0.0001000,
- -0.0005000,
- -0.0010000,
- -0.0050000,
- -0.0100000,
- -0.0500000,
- -0.1000000,
- -0.5000000,
- -1.0000000,
- 0
- };
-
- /* For 2 pass smoothing */
- double smset2[20] = {
- -0.0200000,
- -0.0500000,
- -0.0800000,
- -0.1000000,
- -0.1500000,
- -0.2000000,
- -0.3000000,
- -0.4000000,
- -0.5000000,
- -0.6000000,
- -0.7000000,
- -0.8000000,
- -1.0000000,
- 0
+ double smset[4][20] = {
+ {
+ -0.00000001,
+ -0.00000010,
+ -0.00000100,
+ -0.00001000,
+ -0.00010000,
+ -0.00100000,
+ -0.01000000,
+ -0.10000000,
+ -1.00000000,
+ 0.0
+ },
+ {
+ -0.0000001,
+ -0.0000010,
+ -0.0000100,
+ -0.0001000,
+ -0.0010000,
+ -0.0100000,
+ -0.1000000,
+ -1.0000000,
+ 0.0
+ },
+ {
+ -0.0000001,
+ -0.0000010,
+ -0.0000100,
+ -0.0001000,
+ -0.0010000,
+ -0.0100000,
+ -0.1000000,
+ -1.0000000,
+ 0.0
+ },
+ {
+ -0.0000010,
+ -0.0000100,
+ -0.0001000,
+ -0.0010000,
+ -0.0100000,
+ -0.1000000,
+ -1.0000000,
+ -10.000000,
+ 0.0
+ }
};
-
/* Set of sample points to explore */
int nset[4][20] = {
{
5, 10, 20, 50, 100, 200, 0 /* di = 1 */
},
{
- 25, 100, 400, 2500, 10000, 40000, 0, /* di = 2 */
+ 25, 100, 200, 400, 1000, 2500, 10000, 40000, /* di = 2 */
},
{
- 25, 50, 75, 125, 250, 500, 1000, 2000, 8000, 125000, 0, /* di = 3 */
+ 10, 25, 75, 125, 250, 500, 1000, 2000, 4000, 8000, 16000, 30000, 100000 /* di = 3 */
},
{
- 50, 100, 200, 450, 625, 900, 1800, 3600, 10000, 160000, 1000000, 0, /* di = 4 */
+ 100, 200, 450, 625, 900, 1200, 1800, 3600, 10000, 200000, 500000 /* di = 4 */
}
};
- /* Set of total noise levels to explore */
/* Set of noise levels to explore (average deviation * 4) */
double noiseset[4][20] = {
{
0.0, /* Perfect data */
+ 0.005, /* 0.2 % */
0.01, /* 1.0 % */
0.02, /* 2.0 % */
0.05, /* 5.0 % */
0.10, /* 10.0 % */
- 0.20, /* 20.0 % */
- -1.0,
+// 0.20, /* 20.0 % */
+ -1.0
},
{
0.0, /* Perfect data */
+ 0.005, /* 0.2 % */
0.01, /* 1.0 % */
0.02, /* 2.0 % */
0.05, /* 5.0 % */
0.10, /* 10.0 % */
- 0.20, /* 20.0 % */
- -1.0,
+// 0.20, /* 20.0 % */
+ -1.0
},
{
0.0, /* Perfect data */
+ 0.005, /* 0.2 % */
0.01, /* 1.0 % */
0.02, /* 2.0 % */
0.05, /* 5.0 % */
0.10, /* 10.0 % */
- 0.20, /* 20.0 % */
- -1.0,
+// 0.20, /* 20.0 % */
+ -1.0
},
{
0.0, /* Perfect data */
+ 0.005, /* 0.2 % */
0.01, /* 1.0 % */
0.02, /* 2.0 % */
- 0.03, /* 3.0 % */
0.05, /* 5.0 % */
0.10, /* 10.0 % */
- 0.20, /* 20.0 % */
- -1.0,
+// 0.20, /* 20.0 % */
+ -1.0
},
};
-
printf("Testing underlying smoothness\n");
printf("Profile is '%s'\n",rco->fname);
- if (twopass)
- printf("Two Pass smoothing\n");
- if (extra)
- printf("Extra fitting\n");
-
/* For dimensions */
if (tdi != 0)
sdi = edi = tdi;
@@ -480,7 +483,7 @@ static void do_series_1(refconv* rco, int tdi, int unif, int tntps, int tnlev, i
continue;
}
- printf("No. Sample points %d\n",ntps);
+ printf("\nNo. Sample points %d\n",ntps);
/* For noise levels */
for (j = tnlev; j < 20; j++) {
@@ -500,13 +503,13 @@ static void do_series_1(refconv* rco, int tdi, int unif, int tntps, int tnlev, i
noise = noiseset[di-1][j];
if (noise < 0.0)
break;
- printf("Noise volume %f%%\n",noise * 100.0);
+ printf("\nNoise volume %f%%\n",noise * 100.0);
/* For each channel combination within profile */
for (ix = 0; ; ix++) {
if (set_refconv(rco, ix, di)) {
- DBG(("set_refconv returned nz with ix %f\n",ix));
+ DBG(("set_refconv returned nz with ix %f di %d\n",ix,di));
break;
}
@@ -515,32 +518,29 @@ static void do_series_1(refconv* rco, int tdi, int unif, int tntps, int tnlev, i
/* For smooth factors */
for (k = 0; k < 20; k++) {
- if (twopass)
- smooth = smset2[k];
- else
- smooth = smset[k];
+ smooth = smset[di-1][k];
if (smooth == 0.0) {
DBG(("smooth == 0\n"));
break;
}
- printf("Smooth %9.7f, ",-smooth); fflush(stdout);
+ printf("Smooth %9.8f, ",-smooth); fflush(stdout);
- do_test(rco, &trmse, &tmaxe, &tavge, verb, plot, di, res, ntps, noise, unif, smooth, twopass, extra);
+ do_test(rco, &trmse, &tmaxe, &tavge, verb, plot, di, res, ntps, noise,
+ unif, smooth);
smv[k] = -smooth;
rmse[k] = trmse;
maxe[k] = tmaxe;
printf("maxerr %f, avgerr %f, rmserr %f\n", tmaxe, tavge, trmse);
}
-// bfit = best(k, rmse, smv); /* Best or RMS */
- bfit = best(k, maxe, smv); /* Best of max error */
- printf("Best smoothness = %9.7f, log10 = %4.1f\n",bfit,log10(bfit));
+ bfit = best(k, rmse, smv); /* Best or RMS */
+ printf("Best smoothness = %9.7f, log10 = %4.2f\n",bfit,log10(bfit));
avgbest += log10(bfit);
}
if (ix > 0) {
avgbest /= (double)ix;
- printf("Average best smoothness of %d = %9.7f, log10 = %4.1f\n",ix,pow(10.0,avgbest),avgbest);
+ printf("Average best smoothness of %d = %9.7f, log10 = %4.2f\n",ix,pow(10.0,avgbest),avgbest);
}
}
@@ -550,9 +550,7 @@ static void do_series_1(refconv* rco, int tdi, int unif, int tntps, int tnlev, i
}
/* Verify the current behaviour with test point number and noise volume */
-static void do_series_2(refconv *rco, int di, int unif, int twopass, int extra) {
- int verb = 0;
- int plot = 0;
+static void do_series_2(int verb, int plot, refconv *rco, int di, int unif) {
int res = 0;
int ntps = 0;
double noise = 0.0;
@@ -562,8 +560,8 @@ static void do_series_2(refconv *rco, int di, int unif, int twopass, int extra)
/* Number of trials to do for each dimension */
int trials[4] = {
- 8,
- 8,
+ 12,
+ 10,
8,
5
};
@@ -657,7 +655,7 @@ static void do_series_2(refconv *rco, int di, int unif, int twopass, int extra)
printf("Smooth %9.7f, ",smooth); fflush(stdout);
- do_test(rco, &trmse, &tmaxe, &tavge, verb, plot, di, res, ntps, noise, unif, smooth, twopass, extra);
+ do_test(rco, &trmse, &tmaxe, &tavge, verb, plot, di, res, ntps, noise, unif, smooth);
printf("maxerr %f, avgerr %f, rmserr %f\n", tmaxe, tavge, trmse);
}
@@ -683,10 +681,9 @@ void usage(void) {
fprintf(stderr," -n no Test ""no"" sample points (default 20, 40, 80, 100)\n");
fprintf(stderr," -a amnt Add total level amnt randomness (default 0.0)\n");
fprintf(stderr," -A n Just do the n'th noise level of series\n");
- fprintf(stderr," -2 Use two pass smoothing\n");
- fprintf(stderr," -x Use extra fitting\n");
fprintf(stderr," -s smooth RSPL extra smoothness factor to test (default 1.0)\n");
fprintf(stderr," -g smooth RSPL underlying smoothness factor to test\n");
+ fprintf(stderr," -X ix Select input channel for 1/2D emulation, 0..3\n");
fprintf(stderr," profile.mpp MPP profile to use\n");
exit(1);
}
@@ -701,6 +698,7 @@ int main(int argc, char *argv[]) {
int series = 0;
int unif = 0;
int di = 0; /* Test input dimensions */
+ int ix = 0; /* One off test output channel */
int its = 3; /* Smooth test itterations */
int res = -1;
int ntps = 0;
@@ -708,8 +706,6 @@ int main(int argc, char *argv[]) {
int nlev = 0;
double smooth = 1.0;
double gsmooth = 0.0;
- int twopass = 0;
- int extra = 0;
int smfunc = 0;
double trmse, tavge, tmaxe;
@@ -737,17 +733,17 @@ int main(int argc, char *argv[]) {
if (argv[fa][1] == '?') {
usage();
- } else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') {
+ } else if (argv[fa][1] == 'v') {
verb = 1;
- } else if (argv[fa][1] == 'p' || argv[fa][1] == 'P') {
+ } else if (argv[fa][1] == 'p') {
plot = 1;
- } else if (argv[fa][1] == 'u' || argv[fa][1] == 'U') {
+ } else if (argv[fa][1] == 'u') {
unif = 1;
/* Test series */
- } else if (argv[fa][1] == 'z' || argv[fa][1] == 'Z') {
+ } else if (argv[fa][1] == 'z') {
fa = nfa;
if (na == NULL) usage();
series = atoi(na);
@@ -758,21 +754,21 @@ int main(int argc, char *argv[]) {
smfunc = 1;
/* Dimension */
- } else if (argv[fa][1] == 'd' || argv[fa][1] == 'D') {
+ } else if (argv[fa][1] == 'd') {
fa = nfa;
if (na == NULL) usage();
di = atoi(na);
if (di <= 0 || di > 4) usage();
/* Resolution */
- } else if (argv[fa][1] == 'r' || argv[fa][1] == 'R') {
+ } else if (argv[fa][1] == 'r') {
fa = nfa;
if (na == NULL) usage();
res = atoi(na);
if (res <= 0) usage();
/* Number of sample points */
- } else if (argv[fa][1] == 'n' || argv[fa][1] == 'N') {
+ } else if (argv[fa][1] == 'n') {
fa = nfa;
if (na == NULL) usage();
ntps = atoi(na);
@@ -792,12 +788,6 @@ int main(int argc, char *argv[]) {
nlev = atoi(na);
if (noise < 0) usage();
- } else if (argv[fa][1] == '2') {
- twopass = 1;
-
- } else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') {
- extra = 1;
-
/* Extra smooth factor */
} else if (argv[fa][1] == 's') {
fa = nfa;
@@ -809,9 +799,15 @@ int main(int argc, char *argv[]) {
} else if (argv[fa][1] == 'g') {
fa = nfa;
if (na == NULL) usage();
- smooth = atof(na);
+ gsmooth = atof(na);
if (gsmooth < 0.0) usage();
+ /* Output channel */
+ } else if (argv[fa][1] == 'X') {
+ fa = nfa;
+ if (na == NULL) usage();
+ ix = atoi(na);
+ if (ix < 0 || ix > 3) usage();
} else
usage();
@@ -849,9 +845,9 @@ int main(int argc, char *argv[]) {
if (series > 0) {
if (series == 1)
- do_series_1(&rco, di, unif, ntps, nlev, twopass, extra);
+ do_series_1(verb, plot, &rco, di, unif, ntps, nlev);
else if (series == 2)
- do_series_2(&rco, di, unif, twopass, extra);
+ do_series_2(verb, plot, &rco, di, unif);
else
error("Unknown series %d\n",series);
return 0;
@@ -895,7 +891,9 @@ int main(int argc, char *argv[]) {
} else {
if (verb) {
+ printf("Profile is '%s'\n",rco.fname);
printf("Dimensions %d\n",di);
+ printf("Outpu chan %d\n",ix);
printf("RSPL resolution %d\n",res);
printf("No. Sample points %d (norm %f)\n",ntps, pow((double)ntps, 1.0/di));
printf("Noise volume total %f, == avg. dev. %f\n",noise, 0.25 * noise);
@@ -905,10 +903,16 @@ int main(int argc, char *argv[]) {
printf("Extra smooth %f\n",smooth);
}
+ ix = 0;
+
+ if (set_refconv(&rco, ix, di)) {
+ error("set_refconv returned nz with ix %f\n",ix);
+ }
+
if (gsmooth > 0.0)
- do_test(&rco, &trmse, &tmaxe, &tavge, verb, plot, di, res, ntps, noise, unif, -gsmooth, twopass, extra);
+ do_test(&rco, &trmse, &tmaxe, &tavge, verb, plot, di, res, ntps, noise, unif, -gsmooth);
else
- do_test(&rco, &trmse, &tmaxe, &tavge, verb, plot, di, res, ntps, noise, unif, smooth, twopass, extra);
+ do_test(&rco, &trmse, &tmaxe, &tavge, verb, plot, di, res, ntps, noise, unif, smooth);
printf("Results: maxerr %f, avgerr %f, rmserr %f\n",
tmaxe, tavge, trmse);
@@ -936,9 +940,7 @@ static void do_test(
int ntps, /* Number of sample points */
double noise, /* Sample point noise volume */
int unif, /* nz for uniform noise, else normal */
- double smooth, /* Smoothness to test, +ve for extra, -ve for underlying */
- int twopass, /* Two pass flag */
- int extra /* Extra fit flag */
+ double smooth /* Smoothness to test, +ve for extra, -ve for underlying */
) {
sobol *so; /* Sobol sequence generator */
co *tps = NULL;
@@ -959,12 +961,6 @@ static void do_test(
gres[j] = res;
}
- if (twopass)
- flags |= RSPL_2PASSSMTH;
-
- if (extra)
- flags |= RSPL_EXTRAFIT2;
-
/* Make repeatable by setting random seed before a test set. */
rand32(0x12345678);
@@ -991,6 +987,7 @@ static void do_test(
so->next(so, tps[i].p);
rco->lookup(rco, out, tps[i].p);
+
/* Add randomness to the PCS values */
/* 0.25 * converts total volume to average deviation */
for (f = 0; f < 3; f++) {
@@ -1009,7 +1006,7 @@ static void do_test(
avgdev[2] = 0.25 * noise;
/* Fit to scattered data */
- if (verb) printf("Fitting the scattered data\n");
+ if (verb) printf("Fitting the scattered data, smooth = %f, avgdev = %f\n",smooth,avgdev[0]);
rss->fit_rspl(rss,
flags, /* Non-mon and clip flags */
tps, /* Test points */
@@ -1023,13 +1020,13 @@ static void do_test(
/* Plot out function values */
if (plot) {
int slice;
- printf("Black is target, Red is rspl\n");
+ printf("L*: Black is target, Red is rspl\n");
+ printf("a*: Green is target, Blue is rspl\n");
+ printf("b*: Yellow is target, Purple is rspl\n");
for (slice = 0; slice < (di+1); slice++) {
co tp; /* Test point */
double x[PLOTRES];
- double ya[PLOTRES];
- double yb[PLOTRES];
- double yc[PLOTRES];
+ double yy[6][PLOTRES];
double pp[MXDI], p1[MXDI], p2[MXDI], ss[MXDI];
int n = PLOTRES;
@@ -1060,7 +1057,9 @@ static void do_test(
/* Reference */
rco->lookup(rco, out, pp);
- ya[i] = 0.01 * out[0];
+ yy[0][i] = out[0];
+ yy[2][i] = out[1];
+ yy[4][i] = out[2];
/* RSPL aproximation */
for (j = 0; j < di; j++)
@@ -1068,19 +1067,20 @@ static void do_test(
if (rss->interp(rss, &tp))
tp.v[0] = -0.1;
- yb[i] = tp.v[0];
-
- /* Crude way of setting the scale: */
- yc[i] = 0.0;
- if (i == (n-1))
- yc[0] = 1.0;
+ yy[1][i] = tp.v[0];
+ yy[3][i] = tp.v[1];
+ yy[5][i] = tp.v[2];
+ /* Increment point */
for (j = 0; j < di; j++)
pp[j] += ss[j];
}
/* Plot the result */
- do_plot(x,ya,yb,yc,n);
+ do_plot6(x,yy[0],yy[1],yy[2],yy[3],yy[4],yy[5],n);
+// do_plot(x,yy[0],yy[1],NULL,n);
+// do_plot(x,yy[2],yy[3],NULL,n);
+// do_plot(x,yy[4],yy[5],NULL,n);
}
}
@@ -1091,7 +1091,7 @@ static void do_test(
// so->reset(so);
- /* Fit to scattered data */
+ /* Check fit to scattered data */
if (verb) printf("Fitting the scattered data\n");
for (i = 0; i < 100000; i++) {
// for (i = 0; i < 100; i++) {
@@ -1108,9 +1108,10 @@ static void do_test(
/* RSPL aproximation */
rss->interp(rss, &tp);
-//printf("~1 point %f %f %f -> ref %f %f %f, test %f %f %f\n", tp.p[0], tp.p[1], tp.p[2], out[0], out[1], out[2], tp.v[0], tp.v[1], tp.v[2]);
-
err = icmLabDE(out, tp.v);
+
+//printf("~1 %f: point %f %f %f -> ref %f %f %f, test %f %f %f\n", err, tp.p[0], tp.p[1], tp.p[2], out[0], out[1], out[2], tp.v[0], tp.v[1], tp.v[2]);
+
avge += err;
rmse += err * err;
if (err > maxe)
diff --git a/rspl/smtnd.c b/rspl/smtnd.c
index cdef0a2..f105429 100644
--- a/rspl/smtnd.c
+++ b/rspl/smtnd.c
@@ -26,12 +26,12 @@
#include "rspl.h"
#include "numlib.h"
#include "xicc.h" /* For mpp support */
-#include "plot.h"
#include "rspl_imp.h"
#include "counters.h" /* Counter macros */
+#include "plot.h"
+#include "ui.h"
-/* rspl flags */
-#define MXCHPARAMS 8
+#define MXCHPARAMS 10 /* Input channel curve parameters */
#define PLOTRES 256
@@ -48,14 +48,15 @@ typedef struct {
/* Setup a random function in the given dimensions */
static void setup_func(funcp *p, int di) {
double mn,mx;
+ double ax[MXDI][2];
int i, j;
p->di = di;
- /* Setup random input parameters */
+ /* Setup random input curve parameters */
/* (This is the one that effects smoothness of function the most) */
for (j = 0; j < di; j++) {
- for (mx = 3.0, i = 0; i < MXCHPARAMS; i++, mx *= 0.6) {
+ for (mx = 4.0, i = 0; i < MXCHPARAMS; i++, mx *= 0.6) {
p->ip[j][i] = d_rand(-mx, mx);
}
}
@@ -63,15 +64,32 @@ static void setup_func(funcp *p, int di) {
/* Setup random shape parameters */
for (j = 0; j < di; j++) {
for (i = 0; i < (1 << di); i++) { /* Initially random */
- p->shape[j][i] = d_rand(-1.0, 1.0);
+ p->shape[j][i] = d_rand(-0.1, 0.1);
}
}
- /* Setup the random output parameters */
- mn = 2.0;
- mx = -1.0;
+ /* Setup the random output value parameters */
+
+ /* First some axis dominant values */
+ for (i = 0; i < di; i++) {
+ ax[i][0] = d_rand(0.0, 1.0);
+ ax[i][1] = ax[i][0] + d_rand(-1.0, 1.0);
+ }
+
+ /* Sum them orthogonally and add indepent terms */
+ mn = 5.0;
+ mx = -5.0;
for (i = 0; i < (1 << di); i++) { /* Initially random */
- p->op[i] = d_rand(0.0, 1.0);
+ p->op[i] = 0.0;
+
+ for (j = 0; j < di; j++) {
+ if ((1 << j) & i)
+ p->op[i] += ax[j][1];
+ else
+ p->op[i] += ax[j][0];
+ }
+ p->op[i] += d_rand(-0.3, 0.3);
+
if (p->op[i] < mn)
mn = p->op[i];
if (p->op[i] > mx)
@@ -100,7 +118,7 @@ static double lookup_func(funcp *p, double *v) {
for (m = 0; m < di; m++)
ww[m] = 0.0;
- /* Lookup the shape values */
+ /* Interpolate the shape values */
for (k = 0; k < (1 << di); k++) { /* For each interp vertex */
double vv;
for (vv = 1.0, m = 0; m < di; m++) { /* Compute weighting */
@@ -156,8 +174,8 @@ static void do_test(
double noise, /* Sample point noise volume */
int unif, /* NZ if uniform rather than standard deistribution noise */
double smooth, /* Smoothness to test */
- int twopass, /* Two pass flag */
- int extra /* Extra fit flag */
+ int autosm, /* Use auto smoothing */
+ int seed /* Random seed value offset */
);
/* Compute smoothness of function */
@@ -176,10 +194,10 @@ static double do_stest(
/* Return the optimal smoothness value, based on the */
/* minimum RMS value. */
static double best(int n, double *rmse, double *smv) {
- int i;
+ int i, bi;
rspl *curve;
co *tps = NULL;
- int ns = 500; /* Number of samples */
+ int ns = 2000; /* Number of samples */
datai low,high;
int gres[1];
datai dlow,dhigh;
@@ -214,7 +232,7 @@ static double best(int n, double *rmse, double *smv) {
n, /* Number of test points */
NULL, NULL, gres, /* Low, high, resolution of grid */
NULL, NULL, /* Default data scale */
- -0.00001, /* Underlying smoothing */
+ -0.0007, /* Underlying smoothing */
avgdev, /* Average deviation */
NULL); /* iwidth */
@@ -228,6 +246,45 @@ static double best(int n, double *rmse, double *smv) {
printf("Point %d at %f, should be %f is %f\n",i,log10(smv[i]),rmse[i],tp.v[0]);
}
+#endif
+
+ /* Choose a solution */
+ brmse = 1e38;
+
+ /* Find lowest rms error point */
+ for (i = ns-1; i >= 0; i--) {
+ co tp;
+ double vi;
+
+ vi = i/(ns-1.0);
+ tp.p[0] = log10(smv[0]) + (log10(smv[n-1]) - log10(smv[0])) * vi;
+ curve->interp(curve, &tp);
+
+ if (tp.v[0] < brmse) {
+ blsmv = tp.p[0];
+ brmse = tp.v[0];
+ bi = i;
+ }
+ }
+
+ /* Then increase smoothness until fit error is 1% higher */
+ for (i = bi+1; i < ns; i++) {
+ co tp;
+ double vi;
+
+ vi = i/(ns-1.0);
+ tp.p[0] = log10(smv[0]) + (log10(smv[n-1]) - log10(smv[0])) * vi;
+ curve->interp(curve, &tp);
+
+ if (tp.v[0] >= (1.01 * brmse)) {
+ blsmv = tp.p[0];
+ brmse = tp.v[0];
+ break;
+ }
+ }
+ rv = pow(10.0, blsmv);
+
+#ifdef NEVER
#define TPRES 100
/* Plot the result */
@@ -243,27 +300,11 @@ static double best(int n, double *rmse, double *smv) {
xx[i] = tp.p[0];
yy[i] = tp.v[0];
}
+ printf("Best at %f\n",blsmv);
do_plot(xx,yy,NULL,NULL,TPRES);
}
#endif
- /* Choose a solution */
- brmse = 1e38;
- for (i = 0; i < ns ; i++) {
- co tp;
- double vi;
-
- vi = i/(ns-1.0);
- tp.p[0] = log10(smv[0]) + (log10(smv[n-1]) - log10(smv[0])) * vi;
- curve->interp(curve, &tp);
-
- if (tp.v[0] < brmse) {
- blsmv = tp.p[0];
- brmse = tp.v[0];
- }
- }
-
- rv = pow(10.0, blsmv);
return rv;
}
@@ -274,7 +315,7 @@ static double best(int n, double *rmse, double *smv) {
/* If tdi != 0, just do the given dimension */
/* If tntps != 0, just do the given number of points */
/* If tnlev != 0, just do the given noise level */
-static void do_series_1(int unif, int tdi, int tntps, int tnlev, int twopass, int extra) {
+static void do_series_1(int unif, int tdi, int tntps, int tnlev, int autosm, int seed) {
int verb = 0;
int plot = 0;
int sdi = 1, edi = 4, di;
@@ -288,10 +329,10 @@ static void do_series_1(int unif, int tdi, int tntps, int tnlev, int twopass, in
/* Number of trials to do for each dimension */
int trials[4] = {
+ 12,
+ 10,
8,
- 6,
- 4,
- 3
+ 5
};
/* Resolution of grid for each dimension */
@@ -302,79 +343,52 @@ static void do_series_1(int unif, int tdi, int tntps, int tnlev, int twopass, in
{ 33, 17, 9, 5 }
};
+#ifndef NEVER
/* Set of sample points to explore */
int nset[4][20] = {
{
- 5, 10, 20, 50, 100, 200, /* di = 1 */
+ 5, 10, 20, 50, 100, 200, 400, 800 /* di = 1 */
},
{
- 25, 100, 400, 2500, 10000, 40000, /* di = 2 */
+ 10, 25, 50, 100, 200, 400, 1000, 2500, 10000, 40000 /* di = 2 */
},
{
- 25, 50, 75, 125, 250, 500, 1000, 2000, 8000, 125000, /* di = 3 */
+ 10, 25, 75, 125, 250, 500, 1000, 2000, 4000, 8000, 16000, 30000, 100000 /* di = 3 */
},
{
- 50, 100, 200, 450, 625, 900, 1800, 3600, 10000, 160000, 1000000, /* di = 4 */
+ 100, 200, 450, 625, 900, 1200, 1800, 3600, 10000, 200000, 500000 /* di = 4 */
}
};
-
- /* Set of smoothnesses to explore */
- double smset[4][20] = {
+#else
+ /* Set of sample points to explore */
+ int nset[4][20] = {
{
- -0.0000001,
- -0.0000010,
- -0.0000100,
- -0.0001000,
- -0.0010000,
- -0.0100000,
- -0.1000000,
- -1.0000000,
- 0.0
+ 20, 50, 100, 200 /* di = 1 */
},
{
- -0.0000001,
- -0.0000010,
- -0.0000100,
- -0.0001000,
- -0.0010000,
- -0.0100000,
- -0.1000000,
- -1.0000000,
- 0.0
+ 100, 400, 2500, 10000 /* di = 2 */
},
{
- -0.0000010,
- -0.0000100,
- -0.0001000,
- -0.0010000,
- -0.0100000,
- -0.1000000,
- -1.0000000,
- 0.0
+ 500, 1000, 2000, 8000 /* di = 3 */
},
{
- -0.0000100,
- -0.0001000,
- -0.0010000,
- -0.0100000,
- -0.1000000,
- -1.0000000,
- -10.000000,
- 0.0
+ 50, 900, 1800, 3600, 10000 /* di = 4 */
}
};
-
- /* Set of smoothnesses for twopass smoothing */
- double smset2[4][20] = {
+#endif
+
+ /* Set of smoothnesses to explore */
+ double smset[4][20] = {
{
- -0.0000001,
- -0.0000010,
- -0.0000100,
- -0.0001000,
- -0.0010000,
- -0.0100000,
- -0.1000000,
- -1.0000000,
+ -0.00000001,
+ -0.00000010,
+ -0.00000100,
+ -0.00001000,
+ -0.00010000,
+ -0.00100000,
+ -0.01000000,
+ -0.10000000,
+ -1.00000000,
0.0
},
{
@@ -410,53 +424,79 @@ static void do_series_1(int unif, int tdi, int tntps, int tnlev, int twopass, in
}
};
+#ifndef NEVER
/* Set of noise levels to explore (average deviation * 4) */
double noiseset[4][20] = {
{
0.0, /* Perfect data */
+ 0.005, /* 0.5 % */
0.01, /* 1.0 % */
0.02, /* 2.0 % */
0.05, /* 5.0 % */
0.10, /* 10.0 % */
- 0.20, /* 20.0 % */
- -1.0,
+// 0.20, /* 20.0 % */
+ -1.0
},
{
0.0, /* Perfect data */
+ 0.005, /* 0.5 % */
0.01, /* 1.0 % */
0.02, /* 2.0 % */
0.05, /* 5.0 % */
0.10, /* 10.0 % */
- 0.20, /* 20.0 % */
- -1.0,
+// 0.20, /* 20.0 % */
+ -1.0
},
{
0.0, /* Perfect data */
+ 0.005, /* 0.5 % */
0.01, /* 1.0 % */
0.02, /* 2.0 % */
0.05, /* 5.0 % */
0.10, /* 10.0 % */
- 0.20, /* 20.0 % */
- -1.0,
+// 0.20, /* 20.0 % */
+ -1.0
},
{
0.0, /* Perfect data */
+ 0.005, /* 0.5 % */
0.01, /* 1.0 % */
0.02, /* 2.0 % */
- 0.03, /* 3.0 % */
0.05, /* 5.0 % */
0.10, /* 10.0 % */
- 0.20, /* 20.0 % */
- -1.0,
+// 0.20, /* 20.0 % */
+ -1.0
},
};
+#else
+
+ /* Set of noise levels to explore (average deviation * 4) */
+ double noiseset[4][20] = {
+ {
+ -1.0
+ },
+ {
+ -1.0
+ },
+ {
+ -1.0
+ },
+ {
+ 0.0, /* Perfect data */
+ 0.005, /* 0.2 % */
+ 0.01, /* 1.0 % */
+ 0.02, /* 2.0 % */
+ 0.05, /* 5.0 % */
+ 0.10, /* 10.0 % */
+// 0.20, /* 20.0 % */
+ -1.0
+ },
+ };
+#endif
+
printf("Testing effect of underlying smoothness factors\n");
- if (twopass)
- printf("Two Pass smoothing\n");
- if (extra)
- printf("Extra fitting\n");
/* For dimensions */
if (tdi != 0)
@@ -472,50 +512,47 @@ static void do_series_1(int unif, int tdi, int tntps, int tnlev, int twopass, in
printf("Dimensions %d\n",di);
printf("RSPL resolution %d\n",res);
- /* For noise levels */
- for (j = tnlev; j < 20; j++) { // All noise levels
- double smv[20];
- double rmse[20];
- double bfit;
+ /* For number of sample points */
+ for (i = 0; i < 20; i++) { // All test points
+ int rpts;
+ ntps = nset[di-1][i];
- if (tnlev != 0 && j != tnlev)
+ if (ntps == 0)
break;
- noise = noiseset[di-1][j];
- if (noise < 0.0)
- break;
- printf("\nNoise volume %f%%, average deviation %f%%\n",noise * 100.0, noise * 25.0);
+ if (tntps != 0 && ntps != tntps) /* Skip any not requested */
+ continue;
+
+ /* Make sure at least 100 points are tested */
+ rpts = 1 + 100/ntps;
+ if (rpts > 5)
+ rpts = 5;
+
+ printf("\nNo. Sample points %d, norm %8.2f, total its %d\n",ntps, pow((double)ntps, 1.0/di),its * rpts);
- /* For number of sample points */
- for (i = 0; i < 20; i++) { // All test points
- int rpts;
- ntps = nset[di-1][i];
+ /* For noise levels */
+ for (j = tnlev; j < 20; j++) { // All noise levels
+ double smv[20];
+ double rmse[20];
+ double bfit;
- if (ntps == 0)
+ if (tnlev != 0 && j != tnlev)
break;
- if (tntps != 0 && ntps != tntps) /* Skip any not requested */
- continue;
-
- /* Make sure at least 100 points are tested */
- rpts = 1 + 100/ntps;
- if (rpts > 5)
- rpts = 5;
-
- printf("\nNo. Sample points %d, norm %8.2f, total its %d\n",ntps, pow((double)ntps, 1.0/di),its * rpts);
+ noise = noiseset[di-1][j];
+ if (noise < 0.0)
+ break;
+ printf("\nNoise volume %f%%, average deviation %f%% Log ADev %f\n",noise * 100.0, noise * 25.0, noise > 0.0 ? log10(0.25 * noise) : -9.0);
/* For smooth factors */
for (k = 0; k < 20; k++) { // All smoothing levels
- if (twopass)
- smooth = smset2[di-1][k];
- else
- smooth = smset[di-1][k];
+ smooth = smset[di-1][k];
if (smooth == 0.0)
break;
printf("Underlying smooth %9.7f, ",-smooth); fflush(stdout);
- do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its * rpts, res, ntps, noise, unif,smooth, twopass, extra);
+ do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its * rpts, res, ntps, noise, unif, smooth, autosm, seed);
smv[k] = -smooth;
rmse[k] = trmse;
printf("maxerr %f%%, avgerr %f%%, rmserr %f%%\n",
@@ -523,7 +560,7 @@ static void do_series_1(int unif, int tdi, int tntps, int tnlev, int twopass, in
}
bfit = best(k, rmse, smv);
- printf("Best smoothness = %9.7f, log10 = %4.1f\n",bfit,log10(bfit));
+ printf("Best smoothness = %9.7f, log10 = %4.2f\n",bfit,log10(bfit));
}
}
}
@@ -532,7 +569,7 @@ static void do_series_1(int unif, int tdi, int tntps, int tnlev, int twopass, in
}
/* Explore performance of "optimised" smoothness over test point number and noise volume */
-static void do_series_2(int unif, int twopass, int extra) {
+static void do_series_2(int unif, int autosm, int seed) {
int verb = 0;
int plot = 0;
int di = 0;
@@ -643,7 +680,7 @@ static void do_series_2(int unif, int twopass, int extra) {
double bfit;
noise = noiseset[j];
- printf("Noise volume %f%%, average deviation %f%%\n",noise * 100.0, noise * 25.0);
+ printf("Noise volume %f%%, average deviation %f%% Log ADev %f\n",noise * 100.0, noise * 25.0, noise > 0.0 ? log10(0.25 * noise) : -9.0);
/* For smooth factors */
for (k = 0; k < 5; k++) {
@@ -651,14 +688,14 @@ static void do_series_2(int unif, int twopass, int extra) {
printf("Extra smooth %f, ",smooth); fflush(stdout);
- do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, smooth, twopass, extra);
+ do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, smooth, autosm, seed);
rmse[k] = trmse;
printf("maxerr %f%%, avgerr %f%%, rmserr %f%%\n",
tmaxe * 100.0, tavge * 100.0, trmse * 100.0);
}
bfit = best(5, rmse, smset);
- printf("Best smoothness = %9.7f, log10 = %4.1f\n",bfit,log10(bfit));
+ printf("Best smoothness = %9.7f, log10 = %4.2f\n",bfit,log10(bfit));
}
}
printf("\n");
@@ -676,17 +713,17 @@ void usage(void) {
fprintf(stderr," 1 = Underlying smoothness\n");
fprintf(stderr," 2 = Verify optimised smoothness\n");
fprintf(stderr," -S Compute smoothness factor instead\n");
- fprintf(stderr," -u Use uniformly distributed noise\n");
- fprintf(stderr," -d n Test ""d"" dimension, 1-4 (default 1)\n");
+ fprintf(stderr," -u Use uniformly distributed noise rather than normal\n");
+ fprintf(stderr," -d n Test ""d"" dimension only, 1-4 (default 1)\n");
fprintf(stderr," -t n Test ""n"" random functions (default 1)\n");
fprintf(stderr," -r res Rspl resolution (defaults 129, 65, 33, 17)\n");
fprintf(stderr," -n no Test ""no"" sample points (default 20, 40, 80, 100)\n");
fprintf(stderr," -a amnt Add total randomness to function value (default 0.0)\n");
fprintf(stderr," -A n Just do the n'th noise level of series\n");
- fprintf(stderr," -2 Use two pass smoothing\n");
- fprintf(stderr," -x Use extra fitting\n");
fprintf(stderr," -s smooth RSPL extra smoothness factor to test (default 1.0)\n");
fprintf(stderr," -g smooth RSPL underlying smoothness factor to test\n");
+ fprintf(stderr," -x Use auto smoothing\n");
+ fprintf(stderr," -Z seed Random seed value\n");
exit(1);
}
@@ -695,7 +732,7 @@ int main(int argc, char *argv[]) {
int verb = 0;
int plot = 0;
int series = 0;
- int unif = 0;
+ int unif = 0; /* default normal noise distribution */
int di = 0;
int its = 1;
int res = -1;
@@ -704,10 +741,10 @@ int main(int argc, char *argv[]) {
int nlev = 0;
double smooth = 1.0;
double gsmooth = 0.0;
- int twopass = 0;
- int extra = 0;
+ int autosm = 0;
int smfunc = 0;
double trmse, tavge, tmaxe;
+ int seed = 0;
error_program = "smtnd";
@@ -747,17 +784,17 @@ int main(int argc, char *argv[]) {
if (argv[fa][1] == '?') {
usage();
- } else if (argv[fa][1] == 'v' || argv[fa][1] == 'V') {
+ } else if (argv[fa][1] == 'v' ) {
verb = 1;
- } else if (argv[fa][1] == 'p' || argv[fa][1] == 'P') {
+ } else if (argv[fa][1] == 'p') {
plot = 1;
- } else if (argv[fa][1] == 'u' || argv[fa][1] == 'U') {
+ } else if (argv[fa][1] == 'u') {
unif = 1;
/* Test series */
- } else if (argv[fa][1] == 'z' || argv[fa][1] == 'Z') {
+ } else if (argv[fa][1] == 'z') {
fa = nfa;
if (na == NULL) usage();
series = atoi(na);
@@ -768,29 +805,28 @@ int main(int argc, char *argv[]) {
smfunc = 1;
/* Dimension */
- } else if (argv[fa][1] == 'd' || argv[fa][1] == 'D') {
+ } else if (argv[fa][1] == 'd') {
fa = nfa;
if (na == NULL) usage();
di = atoi(na);
if (di <= 0 || di > 4) usage();
-printf("~1 Got -d %s = %d\n",na,di);
/* Number of tests */
- } else if (argv[fa][1] == 't' || argv[fa][1] == 'T') {
+ } else if (argv[fa][1] == 't') {
fa = nfa;
if (na == NULL) usage();
its = atoi(na);
if (its <= 0) usage();
/* Resolution */
- } else if (argv[fa][1] == 'r' || argv[fa][1] == 'R') {
+ } else if (argv[fa][1] == 'r') {
fa = nfa;
if (na == NULL) usage();
res = atoi(na);
if (res <= 0) usage();
/* Number of sample points */
- } else if (argv[fa][1] == 'n' || argv[fa][1] == 'N') {
+ } else if (argv[fa][1] == 'n') {
fa = nfa;
if (na == NULL) usage();
ntps = atoi(na);
@@ -810,12 +846,6 @@ printf("~1 Got -d %s = %d\n",na,di);
nlev = atoi(na);
if (noise < 0) usage();
- } else if (argv[fa][1] == '2') {
- twopass = 1;
-
- } else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') {
- extra = 1;
-
/* Extra smooth factor */
} else if (argv[fa][1] == 's') {
fa = nfa;
@@ -827,9 +857,18 @@ printf("~1 Got -d %s = %d\n",na,di);
} else if (argv[fa][1] == 'g') {
fa = nfa;
if (na == NULL) usage();
- smooth = atof(na);
+ gsmooth = atof(na);
if (gsmooth < 0.0) usage();
+ } else if (argv[fa][1] == 'x') {
+ autosm = 1;
+
+ /* Random seed offset */
+ } else if (argv[fa][1] == 'Z') {
+ fa = nfa;
+ if (na == NULL) usage();
+ seed = atoi(na);
+
} else
usage();
} else
@@ -838,9 +877,9 @@ printf("~1 Got -d %s = %d\n",na,di);
if (series > 0) {
if (series == 1)
- do_series_1(unif, di, ntps, nlev, twopass, extra);
+ do_series_1(unif, di, ntps, nlev, autosm, seed);
else if (series == 2)
- do_series_2(unif, twopass, extra);
+ do_series_2(unif, autosm, seed);
else
error("Unknown series %d\n",series);
return 0;
@@ -899,9 +938,9 @@ printf("~1 Got -d %s = %d\n",na,di);
}
if (gsmooth > 0.0)
- do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, -gsmooth, twopass, extra);
+ do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, -gsmooth, autosm, seed);
else
- do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, smooth, twopass, extra);
+ do_test(&trmse, &tmaxe, &tavge, verb, plot, di, its, res, ntps, noise, unif, smooth, autosm, seed);
printf("Results: maxerr %f%%, avgerr %f%%, rmserr %f%%\n",
tmaxe * 100.0, tavge * 100.0, trmse * 100.0);
@@ -924,8 +963,8 @@ static void do_test(
double noise, /* Sample point noise volume (total = 4 x average deviation) */
int unif, /* NZ if uniform rather than standard deistribution noise */
double smooth, /* Smoothness to test, +ve for extra, -ve for underlying */
- int twopass, /* Two pass flag */
- int extra /* Extra fit flag */
+ int autosm, /* Use auto smoothing */
+ int seed /* Random seed value offset */
) {
funcp fp; /* Function parameters */
sobol *so; /* Sobol sequence generator */
@@ -937,6 +976,9 @@ static void do_test(
int i, j, it;
int flags = RSPL_NOFLAGS;
+ if (autosm)
+ flags |= RSPL_AUTOSMOOTH;
+
*trmse = 0.0;
*tmaxe = 0.0;
*tavge = 0.0;
@@ -947,20 +989,15 @@ static void do_test(
gres[j] = res;
}
- if (twopass)
- flags |= RSPL_2PASSSMTH;
-
- if (extra)
- flags |= RSPL_EXTRAFIT2;
-
if ((so = new_sobol(di)) == NULL)
error("Creating sobol sequence generator failed");
for (it = 0; it < its; it++) {
double rmse, avge, maxe;
+ double tnoise = 0.0;
/* Make repeatable by setting random seed before a test set. */
- rand32(0x12345678 + 0x1000 * it);
+ rand32(0x12345678 + seed + 0x1000 * it);
/* New function */
setup_func(&fp, di);
@@ -977,16 +1014,24 @@ static void do_test(
if (verb) printf("Generating the sample points\n");
for (i = 0; i < ntps; i++) {
+ double out, n;
so->next(so, tps[i].p);
- tps[i].v[0] = lookup_func(&fp, tps[i].p);
+ out = lookup_func(&fp, tps[i].p);
if (unif)
- tps[i].v[0] += d_rand(-0.5 * noise, 0.5 * noise);
+ n = d_rand(-0.5 * noise, 0.5 * noise);
else
- tps[i].v[0] += noise * 0.25 * 1.2533 * norm_rand();
+ n = noise * 0.25 * 1.2533 * norm_rand();
+
+ tps[i].v[0] = out + n;
+//printf("~1 data %d: %f %f %f -> %f, inc noise %f\n", i, tps[i].p[0], tps[i].p[1], tps[i].p[2], out, tps[i].v[0]);
+
+ tnoise += fabs(n);
}
+ tnoise /= (double) ntps;
+ if (verb) printf("Measured noise average deviation = %f%%\n",tnoise * 100.0);
/* Fit to scattered data */
- if (verb) printf("Fitting the scattered data\n");
+ if (verb) printf("Fitting the scattered data, smooth = %f, avgdev = %f\n",smooth,avgdev != NULL ? avgdev[0] : 0.0);
avgdev[0] = 0.25 * noise;
rss->fit_rspl(rss,
flags, /* Non-mon and clip flags */
@@ -1026,6 +1071,7 @@ static void do_test(
printf("Slice along diagonal\n");
}
+ /* Start point and step increment */
for (j = 0; j < di; j++) {
ss[j] = (p2[j] - p1[j])/n;
pp[j] = p1[j];
diff --git a/rspl/t2d.c b/rspl/t2d.c
index b1167ef..e3f304c 100644
--- a/rspl/t2d.c
+++ b/rspl/t2d.c
@@ -22,6 +22,7 @@
#include "rspl.h"
#include "tiffio.h"
#include "plot.h"
+#include "ui.h"
#ifdef NEVER
@@ -452,8 +453,7 @@ void usage(void) {
fprintf(stderr," -r resx,resy Set grid resolutions (def %d %d)\n",GRES0,GRES1);
fprintf(stderr," -h Test half scale resolution too\n");
fprintf(stderr," -q Test quarter scale resolution too\n");
- fprintf(stderr," -2 Use two pass smoothing\n");
- fprintf(stderr," -x Use extra fitting\n");
+ fprintf(stderr," -x Use auto smoothing\n");
fprintf(stderr," -s Test symetric smoothness (set asymetric -r !)\n");
fprintf(stderr," -S Test spline interpolation\n");
fprintf(stderr," -p plot 3 slices, x = 0.5, y = 0.5, x = y\n");
@@ -473,8 +473,7 @@ int main(int argc, char *argv[]) {
co *test_points = test_points1;
int npoints = sizeof(test_points1)/sizeof(co);
int dospline = 0;
- int twopass = 0;
- int extra = 0;
+ int autosm = 0;
int dosym = 0;
int doplot = 0;
double plotpts[2][2]; /* doplot == 2 start/end points */
@@ -613,11 +612,8 @@ int main(int argc, char *argv[]) {
} else if (argv[fa][1] == 'S') {
dospline = 1;
- } else if (argv[fa][1] == '2') {
- twopass = 1;
-
- } else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') {
- extra = 1;
+ } else if (argv[fa][1] == 'x') {
+ autosm = 1;
} else if (argv[fa][1] == 's') {
dosym = 1;
@@ -632,11 +628,8 @@ int main(int argc, char *argv[]) {
}
- if (twopass)
- flags |= RSPL_2PASSSMTH;
-
- if (extra)
- flags |= RSPL_EXTRAFIT2;
+ if (autosm)
+ flags |= RSPL_AUTOSMOOTH;
if (dosym)
flags |= RSPL_SYMDOMAIN;
diff --git a/rspl/t2ddf.c b/rspl/t2ddf.c
index 92e7e3f..115d5f5 100644
--- a/rspl/t2ddf.c
+++ b/rspl/t2ddf.c
@@ -22,6 +22,7 @@
#include "rspl.h"
#include "tiffio.h"
#include "plot.h"
+#include "ui.h"
#ifdef NEVER
#define INTERP spline_interp
diff --git a/rspl/t3d.c b/rspl/t3d.c
index 1935545..c526669 100644
--- a/rspl/t3d.c
+++ b/rspl/t3d.c
@@ -22,6 +22,7 @@
#include "rspl.h"
#include "tiffio.h"
#include "plot.h"
+#include "ui.h"
#ifdef NEVER
#define INTERP spline_interp
@@ -357,8 +358,7 @@ void usage(void) {
fprintf(stderr," -r resx,resy,resz Set grid resolutions (def %d %d %d)\n",GRES0,GRES1,GRES2);
fprintf(stderr," -h Test half scale resolution too\n");
fprintf(stderr," -q Test quarter scale resolution too\n");
- fprintf(stderr," -2 Use two pass smoothing\n");
- fprintf(stderr," -x Use extra fitting\n");
+ fprintf(stderr," -x Use auto smoothing\n");
fprintf(stderr," -s Test symetric smoothness\n");
fprintf(stderr," -p plot 4 slices, xy = 0.5, yz = 0.5, xz = 0.5, x=y=z\n");
fprintf(stderr," -P x1:y1:z1:x2:y2:z2 plot slice from x1,y1,z1,x2,y2,z2\n");
@@ -377,8 +377,7 @@ int main(int argc, char *argv[]) {
co *test_points = test_points1;
int npoints = sizeof(test_points1)/sizeof(co);
int dosym = 0;
- int twopass = 0;
- int extra = 0;
+ int autosm = 0;
int doplot = 0;
double plotpts[2][3]; /* doplot == 2 start/end points */
int doh = 0; /* half scale */
@@ -482,11 +481,8 @@ int main(int argc, char *argv[]) {
} else if (argv[fa][1] == 's') {
dosym = 1;
- } else if (argv[fa][1] == '2') {
- twopass = 1;
-
- } else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') {
- extra = 1;
+ } else if (argv[fa][1] == 'x') {
+ autosm = 1;
/* smoothing factor */
} else if (argv[fa][1] == 'S') {
@@ -502,11 +498,8 @@ int main(int argc, char *argv[]) {
}
- if (twopass)
- flags |= RSPL_2PASSSMTH;
-
- if (extra)
- flags |= RSPL_EXTRAFIT2;
+ if (autosm)
+ flags |= RSPL_AUTOSMOOTH;
if (dosym)
flags |= RSPL_SYMDOMAIN;
diff --git a/rspl/t3ddf.c b/rspl/t3ddf.c
index b0e7d18..4a43f0d 100644
--- a/rspl/t3ddf.c
+++ b/rspl/t3ddf.c
@@ -22,6 +22,7 @@
#include "rspl.h"
#include "tiffio.h"
#include "plot.h"
+#include "ui.h"
#ifdef NEVER
#define INTERP spline_interp
@@ -103,7 +104,7 @@ void usage(void) {
fprintf(stderr," -r resx,resy,resz Set grid resolutions (def %d %d %d)\n",GRES0,GRES1,GRES2);
fprintf(stderr," -h Test half scale resolution too\n");
fprintf(stderr," -q Test quarter scale resolution too\n");
- fprintf(stderr," -x Use extra fitting\n");
+ fprintf(stderr," -x Use auto smoothing\n");
fprintf(stderr," -s Test symetric smoothness (set asymetric -r !)\n");
fprintf(stderr," -s Test symetric smoothness\n");
fprintf(stderr," -p plot 4 slices, xy = 0.5, yz = 0.5, xz = 0.5, x=y=z\n");
@@ -121,8 +122,7 @@ int main(int argc, char *argv[]) {
co *test_points = test_points1;
int npoints = sizeof(test_points1)/sizeof(co);
double wweight = 1.0;
- int twopass = 0;
- int extra = 0;
+ int autosm = 0;
int dosym = 0;
int doplot = 0;
int doh = 0;
@@ -202,8 +202,8 @@ int main(int argc, char *argv[]) {
} else if (argv[fa][1] == 's' || argv[fa][1] == 'S') {
dosym = 1;
- } else if (argv[fa][1] == 'x' || argv[fa][1] == 'X') {
- extra = 1;
+ } else if (argv[fa][1] == 'x') {
+ autosm = 1;
} else if (argv[fa][1] == 's') {
dosym = 1;
@@ -215,11 +215,8 @@ int main(int argc, char *argv[]) {
}
- if (twopass)
- flags |= RSPL_2PASSSMTH;
-
- if (extra)
- flags |= RSPL_EXTRAFIT2;
+ if (autosm)
+ flags |= RSPL_AUTOSMOOTH;
if (dosym)
flags |= RSPL_SYMDOMAIN;
diff --git a/rspl/tnd.c b/rspl/tnd.c
index 95e4d16..c6f0911 100644
--- a/rspl/tnd.c
+++ b/rspl/tnd.c
@@ -22,6 +22,7 @@
#include "rspl.h"
#include "numlib.h"
#include "tiffio.h"
+#include "plot.h"
#ifdef NEVER
FILE *verbose_out = stdout;
@@ -40,6 +41,10 @@ int verbose_level = 6; /* Current verbosity level */
#undef TEST_SLICE
#undef TEST_RANDOM_POINTS
+#ifdef TEST_SLICE
+# include "ui.h"
+#endif
+
#define MAX_ITS 500
#define IT_TOL 0.0005
#define GRES 17 /* Grid resolution */