summaryrefslogtreecommitdiff
path: root/rspl/scat2.c
diff options
context:
space:
mode:
Diffstat (limited to 'rspl/scat2.c')
-rw-r--r--rspl/scat2.c233
1 files changed, 233 insertions, 0 deletions
diff --git a/rspl/scat2.c b/rspl/scat2.c
new file mode 100644
index 0000000..7579366
--- /dev/null
+++ b/rspl/scat2.c
@@ -0,0 +1,233 @@
+
+/*
+ * Argyll Color Correction System
+ * Multi-dimensional multilevel spline data fitter
+ * mlbs base version.
+ *
+ * Author: Graeme W. Gill
+ * Date: 2000/11/10
+ *
+ * Copyright 1996 - 2000 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.
+ */
+
+/* This file contains the scattered data solution specific code */
+
+/* TTBD:
+ *
+ * mlbs code doesn't work. Results are rubbish.
+ *
+ * Fix bugs ?
+ * merge stest.c into this file.
+ *
+ * Get rid of error() calls - return status instead
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <math.h>
+#include <time.h>
+#if defined(__IBMC__) && defined(_M_IX86)
+#include <float.h>
+#endif
+
+#include "rspl_imp.h"
+#include "numlib.h"
+#include "mlbs.h"
+
+extern void error(char *fmt, ...), warning(char *fmt, ...);
+
+#undef DEBUG
+
+#undef NEVER
+#define ALWAYS
+
+/* Implemented in rspl.c: */
+extern void alloc_grid(rspl *s);
+
+extern int is_mono(rspl *s);
+
+/* ============================================ */
+void set_from_mlbs(void *cbctx, double *out, double *in) {
+ mlbs *p = (mlbs *)cbctx;
+ co tp;
+ int i;
+
+ for (i = 0; i < p->di; i++)
+ tp.p[i] = in[i];
+
+ if (p->lookup(p, &tp))
+ error("Internal, set_from_mlbs failed!");
+
+ for (i = 0; i < p->di; i++)
+ out[i] = tp.v[i];
+}
+
+/* Initialise the regular spline from scattered data */
+/* Return non-zero if non-monotonic */
+static int
+fit_rspl_imp(
+ rspl *s, /* this */
+ int flags, /* Combination of flags */
+ void *d, /* Array holding position and function values of data points */
+ int dtp, /* Flag indicating data type, 0 = (co *), 1 = (cow *) */
+ int dno, /* Number of data points */
+ datai glow, /* Grid low scale - will be expanded to enclose data, NULL = default 0.0 */
+ datai ghigh, /* Grid high scale - will be expanded to enclose data, NULL = default 1.0 */
+ int gres, /* Spline grid resolution */
+ datao vlow, /* Data value low normalize, NULL = default 0.0 */
+ datao vhigh, /* Data value high normalize - NULL = default 1.0 */
+ double smooth /* Smoothing factor, nominal = 1.0 */
+) {
+ int di = s->di, fdi = s->fdi;
+ int i, n, e, f;
+ int rv;
+ int nigc;
+ int bres;
+ mlbs *p;
+
+#if defined(__IBMC__) && defined(_M_IX86)
+ _control87(EM_UNDERFLOW, EM_UNDERFLOW);
+#endif
+ /* set debug level */
+ s->debug = (flags >> 24);
+
+ /* Init other flags */
+ if (flags & RSPL_NONMON) /* Enable elimination of non-monoticities */
+ s->nm = 1;
+ else
+ s->nm = 0;
+
+ if (flags & RSPL_VERBOSE) /* Turn on progress messages to stdout */
+ s->verbose = 1;
+ else
+ s->verbose = 0;
+
+ /* Save smoothing factor */
+ s->smooth = smooth;
+
+ /* Stash the data points away */
+ s->d.no = dno; /* Number of data points */
+
+ /* Allocate the scattered data space */
+ if ((s->d.a = (dpnts *) malloc(sizeof(dpnts) * s->d.no)) == NULL)
+ error("rspl malloc failed - data points");
+
+ if (dtp == 0) { /* Default weight */
+ co *dp = (co *)d;
+
+ /* Copy the list into data points */
+ for (n = 0; n < s->d.no; n++) {
+ for (e = 0; e < s->di; e++)
+ s->d.a[n].p[e] = dp[n].p[e];
+ for (f = 0; f < s->fdi; f++)
+ s->d.a[n].v[f] = dp[n].v[f];
+ s->d.a[n].k = 1.0; /* Assume all data points have same weight */
+ }
+ } else { /* Per data point weight */
+ cow *dp = (cow *)d;
+
+ /* Copy the list into data points */
+ for (n = 0; n < s->d.no; n++) {
+ for (e = 0; e < s->di; e++)
+ s->d.a[n].p[e] = dp[n].p[e];
+ for (f = 0; f < s->fdi; f++)
+ s->d.a[n].v[f] = dp[n].v[f];
+ s->d.a[n].k = dp[n].w; /* Weight specified */
+ }
+ }
+
+ /* Compute target B-Spline resolution */
+ /* Make it worst case half the target rspl resolution */
+ for (bres = 2; (2 * bres) < gres; bres = 2 * bres -1)
+ ;
+
+ /* Create multilevel B-Spline fit */
+ p = new_mlbs(di, fdi, bres, s->d.a, s->d.no, glow, ghigh, smooth);
+ if (p == NULL)
+ error("new_mlbs() failed");
+
+ /* Create rspl grid points by looking up the B-Spline values */
+ rv = s->set_rspl(s, 0, (void *)p, set_from_mlbs, glow, ghigh, gres, vlow, vhigh);
+
+ /* Don't need B-Spline any more */
+ p->del(p);
+
+ return rv;
+}
+
+/* Initialise the regular spline from scattered data */
+/* Return non-zero if non-monotonic */
+int
+fit_rspl(
+ rspl *s, /* this */
+ int flags, /* Combination of flags */
+ co *d, /* Array holding position and function values of data points */
+ int dno, /* Number of data points */
+ datai glow, /* Grid low scale - will be expanded to enclose data, NULL = default 0.0 */
+ datai ghigh, /* Grid high scale - will be expanded to enclose data, NULL = default 1.0 */
+ int gres, /* Spline grid resolution */
+ datao vlow, /* Data value low normalize, NULL = default 0.0 */
+ datao vhigh, /* Data value high normalize - NULL = default 1.0 */
+ double smooth /* Smoothing factor, nominal = 1.0 */
+) {
+ /* Call implementation with (co *) data */
+ return fit_rspl_imp(s, flags, (void *)d, 0, dno, glow, ghigh, gres, vlow, vhigh, smooth);
+}
+
+/* Initialise the regular spline from scattered data with weights */
+/* Return non-zero if non-monotonic */
+int
+fit_rspl_w(
+ rspl *s, /* this */
+ int flags, /* Combination of flags */
+ cow *d, /* Array holding position, function and weight values of data points */
+ int dno, /* Number of data points */
+ datai glow, /* Grid low scale - will be expanded to enclose data, NULL = default 0.0 */
+ datai ghigh, /* Grid high scale - will be expanded to enclose data, NULL = default 1.0 */
+ int gres, /* Spline grid resolution */
+ datao vlow, /* Data value low normalize, NULL = default 0.0 */
+ datao vhigh, /* Data value high normalize - NULL = default 1.0 */
+ double smooth /* Smoothing factor, nominal = 1.0 */
+) {
+ /* Call implementation with (cow *) data */
+ return fit_rspl_imp(s, flags, (void *)d, 1, dno, glow, ghigh, gres, vlow, vhigh, smooth);
+}
+
+/* Init scattered data elements in rspl */
+void
+init_data(rspl *s) {
+ s->d.no = 0;
+ s->d.a = NULL;
+ s->fit_rspl = fit_rspl;
+ s->fit_rspl_w = fit_rspl_w;
+}
+
+/* Free the scattered data allocation */
+void
+free_data(rspl *s) {
+ if (s->d.a != NULL) {
+ free((void *)s->d.a);
+ s->d.a = NULL;
+ }
+}
+
+/* ============================================ */
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+