summaryrefslogtreecommitdiff
path: root/rspl/rspl.c
diff options
context:
space:
mode:
Diffstat (limited to 'rspl/rspl.c')
-rw-r--r--rspl/rspl.c149
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 */