summaryrefslogtreecommitdiff
path: root/rspl/cw3.c
diff options
context:
space:
mode:
Diffstat (limited to 'rspl/cw3.c')
-rw-r--r--rspl/cw3.c325
1 files changed, 325 insertions, 0 deletions
diff --git a/rspl/cw3.c b/rspl/cw3.c
new file mode 100644
index 0000000..03c9b81
--- /dev/null
+++ b/rspl/cw3.c
@@ -0,0 +1,325 @@
+
+/************************************************/
+/* Investigate various curve approximations */
+/************************************************/
+
+/* Discrete regularized spline versions */
+/* Test variable grid spacing in 3D
+
+/* Author: Graeme Gill
+ * Date: 28/11/2013
+ *
+ * Copyright 1995, 2013 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 DIAG
+#undef DIAG2
+#undef GLOB_CHECK
+#define AVGDEV 0.0 /* Average deviation of function data */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <fcntl.h>
+#include <math.h>
+#include "copyright.h"
+#include "aconfig.h"
+#include "counters.h"
+#include "numlib.h"
+#include "plot.h"
+#include "rspl.h"
+
+double lin();
+void usage(void);
+
+#define TRIALS 1 /* Number of random trials */
+#define XRES 100 /* Plotting res */
+
+
+#define PNTS1 20
+#define GRES1 33
+
+/* Function that we're approximating */
+static double func(double val) {
+ double out;
+ int sgn = 0;
+
+ if (val < 0.0) {
+ sgn = 1;
+ val = -val;
+ }
+
+ out = pow(val, 2.4);
+
+ if (sgn)
+ out = -out;
+
+ return out;
+}
+
+/* Grid spacing power */
+#define GRIDPOW 1.8
+
+/* Scattered data sample spacing power */
+#define SAMPPOW 1.3
+
+
+co test_points[PNTS1 * PNTS1 * PNTS1];
+
+double ipos[GRES1];
+double *iposes[3] = { ipos, ipos, ipos };
+
+double powlike(double vv, double pp);
+
+void usage(void) {
+ fprintf(stderr,"Test 3D rspl interpolation variable grid spacing\n");
+ fprintf(stderr,"Author: Graeme W. Gill\n");
+ fprintf(stderr,"usage: c1 [options]\n");
+ fprintf(stderr," -s smooth Use given smoothness (default 1.0)\n");
+ exit(1);
+}
+
+int main(int argc, char *argv[]) {
+ int fa,nfa; /* argument we're looking at */
+ int i, j, e, f, n;
+ double x;
+ double xx1[XRES];
+ double yy1[6][XRES];
+ double xx2[XRES];
+ double yy2[6][XRES];
+ rspl *rss; /* incremental solution version */
+ datai low,high;
+ int gres[MXDI];
+ double smooth = 1.0;
+ double avgdev[MXDO];
+
+ low[0] = 0.0;
+ low[1] = 0.0;
+ low[2] = 0.0;
+ high[0] = 1.0;
+ high[1] = 1.0;
+ high[2] = 1.0;
+ avgdev[0] = AVGDEV;
+ avgdev[1] = AVGDEV;
+ avgdev[2] = AVGDEV;
+
+ error_program = "c1";
+ check_if_not_interactive();
+
+ /* Process the arguments */
+ for(fa = 1;fa < argc;fa++) {
+ nfa = fa; /* skip to nfa if next argument is used */
+ if (argv[fa][0] == '-') { /* Look for any flags */
+ char *na = NULL; /* next argument after flag, null if none */
+
+ if (argv[fa][2] != '\000')
+ na = &argv[fa][2]; /* next is directly after flag */
+ else {
+ if ((fa+1) < argc) {
+ if (argv[fa+1][0] != '-') {
+ nfa = fa + 1;
+ na = argv[nfa]; /* next is seperate non-flag argument */
+ }
+ }
+ }
+
+ if (argv[fa][1] == '?')
+ usage();
+
+ /* smoothness */
+ else if (argv[fa][1] == 's' || argv[fa][1] == 'S') {
+ fa = nfa;
+ if (na == NULL) usage();
+ smooth = atof(na);
+ }
+ else
+ usage();
+ } else
+ break;
+ }
+
+ for (n = 0; n < TRIALS; n++) {
+ double lrand = 0.0; /* Amount of level randomness */
+ int di = 3;
+ int pnts = PNTS1;
+ int tpnts = PNTS1 * PNTS1 * PNTS1;
+ int fres;
+
+ fres = GRES1;
+ gres[0] = fres;
+ gres[1] = fres;
+ gres[2] = fres;
+
+ /* Create the object */
+ rss = new_rspl(RSPL_NOFLAGS,
+ 3, /* di */
+ 3); /* fdi */
+
+ /* For each grid spacing test */
+ for (j = 0; j < 2; j++) {
+ DCOUNT(gc, MXDIDO, di, 0, 0, pnts);
+
+ /* Set the scattered data points (actually grid) */
+ DC_INIT(gc);
+ for (i = 0; !DC_DONE(gc); i++) {
+
+ for (f = 0; f < 3; f++) {
+ double val, out;
+
+ val = gc[f]/(pnts - 1.0);
+
+ if (j == 0) { /* Linear */
+// out = func(val + d_rand(-0.001, 0.001));
+ out = func(val);
+ test_points[i].p[f] = val;
+ test_points[i].v[f] = out;
+
+ } else { /* power 2.0 grid spacing */
+ val = pow(val, SAMPPOW); /* sample point spacing */
+
+// out = func(val + d_rand(-0.001, 0.001));
+ out = func(val);
+ val = powlike(val, 1.0/GRIDPOW); /* To grid spacing */
+ test_points[i].p[f] = val;
+ test_points[i].v[f] = out;
+ }
+ }
+ DC_INC(gc);
+ }
+
+ /* Set grid position info. This is the same for each input axis */
+ for (i = 0; i < gres[0]; i++) {
+ double val, out;
+
+ val = i/(gres[0] - 1.0);
+ if (j == 1)
+ val = powlike(val, 1.0/GRIDPOW); /* To grid spacing */
+ ipos[i] = val;
+ }
+
+ printf("Fitting rspl....\n");
+ rss->fit_rspl(rss,
+ 0,
+ test_points, /* Test points */
+ tpnts, /* Number of test points */
+ low, high, gres, /* Low, high, resolution of grid */
+ NULL, NULL, /* Default data scale */
+ smooth, /* Smoothing */
+ avgdev, /* Average deviation */
+ j == 0 ? NULL : iposes); /* iwidth */
+
+ /* Show the result full scale of one axis and corresponding output */
+ for (i = 0; i < XRES; i++) {
+ co tp; /* Test point */
+ x = i/(double)(XRES-1);
+ xx1[i] = x;
+ yy1[0][i] = func(x);
+ tp.p[0] = x;
+ tp.p[1] = tp.p[2] = 0.0;
+ if (j == 1)
+ tp.p[0] = powlike(tp.p[0], 1.0/GRIDPOW); /* Grid spacing input */
+ rss->interp(rss, &tp);
+ yy1[1+j][i] = tp.v[0];
+ if (yy1[1+j][i] < -0.2)
+ yy1[1+j][i] = -0.2;
+ else if (yy1[1+j][i] > 1.2)
+ yy1[1+j][i] = 1.2;
+ }
+
+ /* Show the result magnified */
+ for (i = 0; i < XRES; i++) {
+ co tp; /* Test point */
+ x = i/(double)(XRES-1);
+ x *= 0.1;
+ xx2[i] = x;
+ yy2[0][i] = func(x);
+ tp.p[0] = x;
+ tp.p[1] = tp.p[2] = 0.0;
+ if (j == 1)
+ tp.p[0] = powlike(tp.p[0], 1.0/GRIDPOW); /* Grid spacing input */
+ rss->interp(rss, &tp);
+ yy2[1+j][i] = tp.v[0];
+ if (yy2[1+j][i] < -0.2)
+ yy2[1+j][i] = -0.2;
+ else if (yy2[1+j][i] > 1.2)
+ yy2[1+j][i] = 1.2;
+ }
+ }
+
+ printf("Full scale: Black = lin, Red = even, Green = pow2\n");
+ do_plot6(xx1,yy1[0],yy1[1],yy1[2],NULL,NULL,NULL,XRES);
+
+ printf("Magnified: Black = lin, Red = even, Green = pow2\n");
+ do_plot6(xx2,yy2[0],yy2[1],yy2[2],NULL,NULL,NULL,XRES);
+ }
+ return 0;
+}
+
+/* A power-like function, based on Graphics Gems adjustment curve. */
+/* Avoids "toe" problem of pure power. */
+/* Adjusted so that "power" 2 and 0.5 agree with real power at 0.5 */
+
+double powlike(double vv, double pp) {
+ double tt, g;
+
+ if (pp >= 1.0) {
+ g = 2.0 * (pp - 1.0);
+ vv = vv/(g - g * vv + 1.0);
+ } else {
+ g = 2.0 - 2.0/pp;
+ vv = (vv - g * vv)/(1.0 - g * vv);
+ }
+
+ return vv;
+}
+
+
+/******************************************************************/
+/* Error/debug output routines */
+/******************************************************************/
+
+/* Next u function done with optimization */
+
+/* Structure to hold data for optimization function */
+struct _edatas {
+ rspl *rss;
+ int j;
+ }; typedef struct _edatas edatas;
+
+#ifdef GLOB_CHECK
+/* Overall Global optimization method */
+/* Definition of the optimization function handed to powell() */
+double efunc2(void *edata, double p[])
+ {
+ int j;
+ double rv;
+ rspl *rss = (rspl *)edata;
+ for (j = 0; j < rss->nig; j++) /* Ugg */
+ rss->u[j].v = p[j];
+ rv = rss->efactor(rss);
+#ifdef DIAG2
+ /* printf("%c%e",cr_char,rv); */
+ printf("%e\n",rv);
+#endif
+ return rv;
+ }
+
+solveu(rss)
+rspl *rss;
+ {
+ int j;
+ double *cp;
+ double *s;
+
+ cp = dvector(0,rss->nig);
+ s = dvector(0,rss->nig);
+ for (j = 0; j < rss->nig; j++) /* Ugg */
+ {
+ cp[j] = rss->u[j].v;
+ s[j] = 0.1;
+ }
+ powell(rss->nig,cp,s,1e-7,1000,efunc2,(void *)rss);
+ }
+#endif /* GLOB_CHECK */
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 17:06:47 +0000