diff options
Diffstat (limited to 'rspl/trnd.c')
| -rw-r--r-- | rspl/trnd.c | 275 |
1 files changed, 275 insertions, 0 deletions
diff --git a/rspl/trnd.c b/rspl/trnd.c new file mode 100644 index 0000000..2676182 --- /dev/null +++ b/rspl/trnd.c @@ -0,0 +1,275 @@ + +/************************************************/ +/* Test RSPL reverse lookup in 3/4D */ +/************************************************/ + +/* Author: Graeme Gill + * Date: 31/10/96 + * Derived from tnd.c + * Copyright 1999 - 2000 Graeme W. Gill + * + * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :- + * see the License.txt file for licencing details. + */ + + +#undef DEBUG +#undef DETAILED + +#include <stdio.h> +#include <fcntl.h> +#include <math.h> +#include "rspl.h" +#include "numlib.h" + +#ifdef NEVER +FILE *verbose_out = stdout; +int verbose_level = 6; /* Current verbosity level */ + /* 0 = none */ + /* !0 = diagnostics */ +#endif /* NEVER */ + +#define GRES 10 /* Grid resolution */ +#define DI 4 /* Dimensions in */ +#define FDI 3 /* Function (out) Dimensions */ +#define DOLIMIT +#define LIMITV 1.50 + +/* Fwd function approximated by rspl */ +void func( +void *cbctx, +double *out, +double *in) { + double tt[4]; + +#ifdef NEVER +printf(" Got input %f %f\n",in[0],in[1]); + if (in[1] < 0.5 && in[0] < 0.5) { /* 0,0 */ + out[0] = 0.1; + out[1] = 0.5; + } + if (in[1] < 0.5 && in[0] > 0.5) { /* 0,1 */ + out[0] = 0.5; + out[1] = 0.0; + } + if (in[1] > 0.5 && in[0] < 0.5) { /* 0,1 */ + out[0] = 0.9; + out[1] = 0.8; + } + if (in[1] > 0.5 && in[0] > 0.5) { /* 0,1 */ + out[0] = 0.9; + out[1] = 0.2; + } +#endif /* NEVER */ + +#if DI >= 3 + tt[0] = 0.7 * in[0] + 0.2 * in[1] + 0.1 * in[2]; + tt[0] = tt[0] * tt[0]; + tt[1] = 0.2 * in[0] + 0.8 * in[1] - 0.1 * in[2]; + tt[1] = tt[1] * tt[1]; + tt[2] = 0.3 * in[0] - 0.2 * in[1] + 0.9 * in[2]; + tt[2] = tt[2] * tt[2]; +#endif + +#if DI == 4 + tt[0] *= in[3]; + tt[1] *= in[3]; + tt[2] *= in[3]; +#endif + + tt[3] = 0.3 * in[0] + 0.4 * in[1] + 0.3 * in[2]; + +#if FDI > 0 + out[0] = tt[0]; +#endif +#if FDI > 1 + out[1] = tt[1]; +#endif +#if FDI > 2 + out[2] = tt[2]; +#endif +#if FDI > 3 + out[3] = tt[3]; +#endif +} + +/* Simplex ink limit function */ +double limitf( +void *lcntx, +double *in +) { + int i; + double ov; + + for (ov = 0.0, i = 0; i < DI; i++) { + ov += in[i]; + } + return ov; +} + +int +main( +int argc, +char *argv[] +) { + int e; + rspl *rss; /* Multi-resolution regularized spline structure */ + int gres[MXDI]; + + for (e = 0; e < DI; e++) + gres[e] = GRES; + + printf("Started test\n"); + /* Create the object */ + rss = new_rspl(RSPL_NOFLAGS, DI, FDI); + + printf("Rspl allocated\n"); + rss->set_rspl(rss, 0, (void *)NULL, func, NULL, NULL, gres, NULL, NULL); +// rss->set_rspl(rss, RSPL_SET_APXLS, (void *)NULL, func, NULL, NULL, gres, NULL, NULL); + + printf("Rspl set\n"); + + { +#define NIP 10 /* Number of solutions allowed */ + int i, r, cl; + double v[MXDO]; /* Target output value */ + co tp[NIP], chp; /* Test point, check point */ + double cvec[4]; /* Text clip vector */ + int auxm[4]; /* Auxiliary target value valid flag */ + double lmin[4], lmax[4]; /* Locus min/max values */ + + /* Output value being looked for */ +/* + tp[0].v[0] = v[0] = 1.5; + tp[0].v[1] = v[1] = 0.9; + tp[0].v[2] = v[2] = 1.2; + tp[0].v[3] = v[3] = 0.0; +*/ + tp[0].v[0] = v[0] = 0.3; + tp[0].v[1] = v[1] = 0.4; + tp[0].v[2] = v[2] = 0.26; + tp[0].v[3] = v[3] = 0.0; + + /* Set auxiliary target */ + auxm[0] = 0; + auxm[1] = 0; + auxm[2] = 0; + auxm[3] = 1; + tp[0].p[0] = 0; + tp[0].p[1] = 0; + tp[0].p[2] = 0.0; + tp[0].p[3] = 0.87; + + for (i = 1; i < NIP; i++) { /* Make sure we can see changes */ + tp[i].p[0] = -1.0; + tp[i].p[1] = -1.0; + tp[i].p[2] = -1.0; + tp[i].p[3] = -1.0; + } + + /* Clip center */ + cvec[0] = 0.1 - tp[0].v[0]; + cvec[1] = 0.1 - tp[0].v[1]; + cvec[2] = 0.1 - tp[0].v[2]; + cvec[3] = 0.1 - tp[0].v[3]; + +#ifdef DOLIMIT + /* Setup limit */ + rss->rev_set_limit(rss, + limitf, /* limit function */ + NULL, /* Function context */ + LIMITV /* limit maximum value */ + ); +#endif /* DOLIMIT */ + +#if DI > FDI + /* Check out the locus size */ + if ((r = rss->rev_locus(rss, + auxm, /* auxm Auxiliary mask flags */ + tp, /* Input and auxiliary values */ + lmin, /* Locus min/max return values */ + lmax + )) > 0) { + + printf("Total of %d Results\n",r); + for (i = 0; i < FDI; i++) { + if (auxm[i] == 0) + continue; + printf("Auxiliary %d min = %f, max = %f\n",i, lmin[i], lmax[i]); + } + } else { + printf("Failed to find gamut range\n"); + } + printf("\n\n"); +#endif + + /* Do reverse interpolation ~~~1 */ + if ((r = rss->rev_interp(rss, + 0 /* RSPL_EXACTAUX */, /* Hint flags */ + NIP, /* Number of solutions allowed */ + auxm, /* auxm Auxiliary mask flags */ + cvec, /* cvec Clip vector direction & length */ + tp) /* Input and output values */ + ) > 0) { + + printf("Total of %d Results\n", r &= RSPL_NOSOLNS); + + printf("Target output %f, %f, %f, %f\n", v[0], v[1], v[2], v[3]); + + cl = r & RSPL_DIDCLIP; /* clipped flag */ + r &= RSPL_NOSOLNS; /* Get number of solutions */ + + /* Check test result */ + for (i = 0; i < r; i++) { + chp.p[0] = tp[i].p[0]; + chp.p[1] = tp[i].p[1]; + chp.p[2] = tp[i].p[2]; + chp.p[3] = tp[i].p[3]; + chp.v[0] = -1.0; + chp.v[1] = -1.0; + chp.v[2] = -1.0; + chp.v[3] = -1.0; + printf("Result %d = inp: %f, %f, %f, %f\n",i, tp[i].p[0],tp[i].p[1],tp[i].p[2],tp[i].p[3]); + printf(" out: %f, %f, %f, %f\n",tp[i].v[0],tp[i].v[1],tp[i].v[2],tp[i].v[3]); + + if (rss->interp(rss, &chp)) + printf("Fwd check %d failed!\n",i); + else { + int p; + double er = 0.0; + for (p = 0; p < FDI; p++) + er += (v[p] - chp.v[p]) * (v[p] - chp.v[p]); + er = sqrt(er); + printf("Fwd check error %d = %f\n",i,er); + + if (cl != 0) { + /* Check if clipped result us reasonable */ + /* ~~~ */ + } + } + { + printf("Output limit = %f\n",limitf(NULL, tp[i].p)); + } + } + } else + printf("Rev lookup result returned none\n"); + } + rss->del(rss); + return 0; +} + + +#ifdef NEVER +/* Standard interface for powell function */ +/* return err on sucess, -1.0 on failure */ +/* Result will be in cp */ +double powell( +int di, /* Dimensionality */ +double cp[], /* Initial starting point */ +double s[], /* Size of initial search area */ +double ftol, /* Tollerance of error change to stop on */ +int maxit, /* Maximum iterations allowed */ +double (*func)(void *fdata, double tp[]), /* Error function to evaluate */ +void *fdata /* Opaque data needed by function */ +) { +#endif |