New release 1.7.0

1 files changed, 149 insertions, 0 deletions

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 */