From 22f703cab05b7cd368f4de9e03991b7664dc5022 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Mon, 1 Sep 2014 13:56:46 +0200 Subject: Initial import of argyll version 1.5.1-8 --- rspl/spline.c | 352 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 352 insertions(+) create mode 100644 rspl/spline.c (limited to 'rspl/spline.c') diff --git a/rspl/spline.c b/rspl/spline.c new file mode 100644 index 0000000..dabeb0f --- /dev/null +++ b/rspl/spline.c @@ -0,0 +1,352 @@ + +/* + * Argyll Color Correction System + * Multi-dimensional regularized spline data structure + * + * Spline forward interpolation support. + * + * Author: Graeme W. Gill + * Date: 12/10/98 + * + * Copyright 1998, Graeme W. Gill + * All rights reserved. + * + * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :- + * see the License.txt file for licencing details. + */ + +/* TTBD: + Get rid of error() calls - return status instead + */ + +#include +#include +#include +#include +#include +//#if defined(__IBMC__) && defined(_M_IX86) +//#include +//#endif + +#include "rspl_imp.h" +#include "numlib.h" + +int spline_interp_rspl(rspl *ss, co *cp); + +#undef DEBUG + +#undef NEVER +#define ALWAYS + +/* Convention is to use: + i to index grid points u.a + n to index data points d.a + e to index position dimension di + f to index output function dimension fdi + j misc and cube corners + k misc + */ + +/* ====================================================== */ + +/* Init spline elements in rspl */ +void +init_spline(rspl *s) { + s->spline.nm = 0; + s->spline.spline = 0; + s->spline.magic = NULL; + + s->spline_interp = spline_interp_rspl; +} + +/* Free up the spline interpolation info */ +void free_spline( +rspl *s /* Pointer to rspl grid */ +) { + if (s->spline.magic != NULL) { + free(s->spline.magic); + } + s->spline.nm = 0; + s->spline.spline = 0; +} + +/* ====================================================== */ +/* Setup functions first: */ + +/* Hermite spline, magic matrix */ +/* Indexes are: param powers 0, 1, 2, 3; Offset from base vertex 0,1; Dimension mask 0,1 */ +static double hmagic[4][2][2] = { + { { 1.0, 0.0}, { 0.0, 0.0} }, + { { 0.0, 1.0}, { 0.0, 0.0} }, + { {-3.0, -2.0}, { 3.0, -1.0} }, + { { 2.0, 1.0}, {-2.0, 1.0} } +}; + +/* Allocate and initialize tangency information for each grid point */ +static void make_tang( +rspl *s /* Pointer to rspl grid */ +) { + int i,p,j; + int di = s->di; + int fdi = s->fdi; + int nig = s->g.no; + float *tp; /* Pointer to tangent values */ + int nim, mix; /* Number in magic, magic index */ + float *tang_alloc, *tang; /* Tangency info */ + float *gt; /* Working grid point */ + +//printf("~~make_tang called\n"); + /* Organized as: tang[[grid]][di combs.][fdi] */ + /* Allocate space for tangency info */ + if ((tang_alloc = (float *) malloc(sizeof(float) * nig * (((1 << di) * fdi)+G_XTRA))) == NULL) + error("rspl malloc failed - tangecy points"); + tang = tang_alloc + G_XTRA; /* Offset for flags and non-mono error */ + + /* For all grid points */ + for (tp = tang, gt = s->g.a, i = 0; i < nig; i++, gt += s->g.pss, tp += G_XTRA) { + int ee; +/* printf("\n~~ grid point %d\n",i); */ + + /* Look at surrounding grid points in combinations of +- 1 all dimensions */ + for (ee = 0; ee < (1 << di); ee++) { + double av[MXRO]; /* average */ + int nia = 0; /* Number in average */ + int f, ec; + +/* printf("Dim combo %d\n",ee); */ + /* special case - base value */ + if (ee == 0) { + *((int *)(tp-2)) = *((int *)(gt-2)); /* Copy flags */ + tp[-1] = gt[-1]; /* Copy ink limit function value */ + for (f = 0; f < fdi; f++) { + *tp++ = gt[f]; +/* printf("Tang value out %d = %f\n",f,tp[-1]); */ + } + continue; + } + for (f = 0; f < fdi; f++) + av[f] = 0.0; /* Init average */ + /* For all surroundin grid points in this combination */ + for (ec = 0; ec < (1 << di); ec++) { + int xo, io, sgn, e, ex; +/* printf("~~checking out surrounding combo %d\n",ec); */ + if (ec & ~ee) { +/* printf("~~being skipped\n"); */ + continue; /* Skip invalid combo */ + } + xo = io = 0; /* Grid float offset */ + sgn = 1; /* Sign */ + ex = 0; /* Flag - No extrapolation */ + for (e = 0; e < di; e++) { /* For each dimension */ +/* printf("~~checking dimension %d\n",e); */ + if (!(ee & (1 << e))) { +/* printf("~~dimension not active\n"); */ + continue; /* Dimension is not active */ + } + if (ec & (1 << e)) { + /* If + dimension is valid */ + if (((G_FL(gt,e) & 3) > 0) || (G_FL(gt,e) & 0x4)) { + int to = s->g.fci[e]; /* +1 in dimension */ + io += to; /* real/pivot point */ + xo += to; /* reflected point */ + } else { + ex = 1; /* Use extrapolation */ + xo -= s->g.fci[e]; /* -1 in dimension */ + } + } else { + sgn = -sgn; /* Reverse sign */ + /* If - dimension is valid */ + if (((G_FL(gt,e) & 3) > 0) || !(G_FL(gt,e) & 0x4)) { + int to = -s->g.fci[e]; /* -1 in dimension */ + io += to; /* real/pivot point */ + xo += to; /* reflected point */ + } else { + ex = 1; /* Use extrapolation */ + xo += s->g.fci[e]; /* +1 in dimension */ + } + } + } + /* Add surrounding grid points value into the average */ + if (!ex) { + for (f = 0; f < fdi; f++) + av[f] += (double)sgn * gt[io + f]; + } else { /* Extrapolate point beyond edge */ + /* Use an extrapolation that tries to maintain curvature */ + for (f = 0; f < fdi; f++) { + double v0,v1,v2; + v0 = gt[io + f]; /* Pivot point */ + v1 = gt[xo + f]; /* Reflection of target in pivot */ + v2 = gt[2 * xo - io + f]; /* Reflection +2 */ + av[f] += (double)sgn * (3.0 * (v0 - v1) + v2); + } + } + nia++; + } + for (f = 0; f < fdi; f++) { + *tp++ = (float)(av[f]/(double)nia); +/* printf("Tang value out %d = %f, average of %d\n",f,tp[-1],nia); */ + } + } /* Next dimension combination */ + } /* Next grid point */ + + /* Create a full sized hermite magic matrix */ + /* Organized as: magic[4^di][2^di][2^di] */ + /* = [param power combos][cube vertex index][di combos], */ + /* but then only store non-zero weight values. */ + for (i = 0, nim = 1; i < di; nim *= 10, i++); /* Number of entries needed */ + if (s->spline.magic == NULL) { /* Allocate space for magic matrix info */ + if ((s->spline.magic = (magic_data *) malloc(sizeof(magic_data) * nim)) == NULL) + error("rspl malloc failed - hermite magic matrix data"); + } + + mix = 0; + for (p = 0; p < (1 << (2 * di)); p++) { /* For all combinations of parameter powers */ + for (i = 0; i < (1 << di); i++) { /* For all corners of cube */ + for (j = 0; j < (1 << di); j++) { /* For all dimension combinations */ + int ii; + double wgt = 1.0; + for (ii = 0; ii < di; ii++) { + wgt *= hmagic[3&(p>>(2*ii))][1&(i>>ii)][1&(j>>ii)]; + } + if (wgt != 0.0) { /* record non-zero weight value */ + s->spline.magic[mix].p = p; + s->spline.magic[mix].i = i; + s->spline.magic[mix].j = fdi * j; /* Pre-scale */ + s->spline.magic[mix].wgt = (float)wgt; + mix++; + } + } + } + } + /* mix should == nim! */ + s->spline.nm = nim; + + /* Free basic grid info, and substitute tangency enhanced version */ + /* ~~~~!! need to free any other structures in rspl that depend on */ + /* ~~~~!! g.pss size, ie. rev stuff ??? */ + if (s->g.alloc != NULL) + free((void *)s->g.alloc); + + s->g.alloc = tang_alloc; + s->g.a = tang; + + /* Adjust index tables */ + s->g.pss = (1 << di) * fdi + G_XTRA; + for (i = 0; i < di; i++) + s->g.fci[i] = s->g.ci[i] * s->g.pss; /* In floats */ + for (i = 0; i < (1 << di); i++) + s->g.fhi[i] = s->g.hi[i] * s->g.pss; /* In floats */ + + s->spline.spline = 1; + +//printf("~~make_tang finished\n"); +} + +/* Do a Hermite spline smooth interpolation based on the finest grid */ +/* (To do this more accurately, the data point interpolation within */ +/* the grid itteration should be of the same order. This increases */ +/* itteration complexity quite a bit, so we won't bother for the moment.) */ +/* This code is not optimised for speed. */ +/* Return 0 if OK, 1 if input was clipped to grid */ +int spline_interp_rspl( +rspl *s, +co *cp /* Input value and returned function value */ +) { + int e,f,p,i; + int di = s->di; + int fdi = s->fdi; + double ppw[MXRI][4]; /* Parameter powers of 0, 1, 2, 3 */ + float *ga[POW2MXRI]; /* Pointers to grid cubes data in tang[] */ + magic_data *tp; /* Pointer to items in magic matrix */ + int rv = 0; + +/* printf("~~smooth interp called\n"); */ + + /* This is a restricted size function */ + if (di > MXRI) + error("rspl: spline can't handle di = %d",di); + if (fdi > MXRO) + error("rspl: spline can't handle fdi = %d",fdi); + + if (s->spline.spline == 0) /* Compute tangent info if it doesn't exist */ + make_tang(s); + + /* Locate grid base point, and position with base cube */ + ga[0] = s->g.a; /* Base pointer of cube */ + for (e = 0; e < di; e++) { + double t, pe; + int mi, gres_1 = s->g.res[e]-1; + + pe = cp->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; + ga[0] += s->g.fci[e] * mi; /* Add offset in dimen */ + t = t - (double)mi;; /* sub-cube offset = parameter in dimension e */ + ppw[e][0] = 1.0; /* Powers of parameter */ + ppw[e][1] = t; + ppw[e][2] = t * t; + ppw[e][3] = t * t * t; + } + + /* Compute indexes into cube corners in tangent array */ + for (i = 1; i < (1 << di); i++) + ga[i] = ga[0] + s->g.fhi[i]; + + /* Now compute the output values */ + for (f = 0; f < fdi; f++) /* Zero output value sums */ + cp->v[f] = 0.0; + + /* For all non-zero combinations of parameter powers */ + { + double ppc = -1000.0; /* Parameter power combination */ + for (tp = s->spline.magic, p = -1; tp < &s->spline.magic[s->spline.nm]; tp++) { + double wgt; /* Magic matrix weight */ + float *gp; /* Pointer to vertex data */ + + if (p != tp->p) { /* Param power needs re-calculating */ + int pp; + p = tp->p; + for (ppc = 1.0, pp = 0; pp < di; pp++) + ppc *= ppw[pp][3&(p>>(2*pp))]; /* comb. of param powers value */ + } + + wgt = tp->wgt * ppc; /* matrix times parameter */ + gp = ga[tp->i] + tp->j; /* Point to base of vertex data */ + for (f = 0; f < fdi; f++) /* For all output values */ + cp->v[f] += wgt * gp[f]; + } + } +/* printf("~~smooth interp finished\n"); */ + return rv; +} + + + + + + + + + + + + + + + + + + + -- cgit v1.2.3