summaryrefslogtreecommitdiff
path: root/rspl
diff options
context:
space:
mode:
authorJörg Frings-Fürst <debian@jff-webhosting.net>2014-09-07 13:29:54 +0200
committerJörg Frings-Fürst <debian@jff-webhosting.net>2014-09-07 13:29:54 +0200
commita879b4e708b3e46c9697ba6581687eeb5b02a320 (patch)
treef85acc5cb1a7c0b03a050c879c20cab3deb9f40e /rspl
parent556dffcdad42b938bc196819aa463247de709765 (diff)
parentc07d0c2d2f6f7b0eb6e92cc6204bf05037957e82 (diff)
state from 2014-09-07 13:30 MESZ
Diffstat (limited to 'rspl')
-rw-r--r--rspl/Jamfile8
-rw-r--r--rspl/Makefile.am15
-rw-r--r--rspl/afiles2
-rw-r--r--rspl/c1.c3
-rw-r--r--rspl/cw1.c308
-rw-r--r--rspl/cw3.c325
-rw-r--r--rspl/rev.c34
-rw-r--r--rspl/rspl.c2
-rw-r--r--rspl/scat.c270
9 files changed, 752 insertions, 215 deletions
diff --git a/rspl/Jamfile b/rspl/Jamfile
index a3f2e8c..dbca6d5 100644
--- a/rspl/Jamfile
+++ b/rspl/Jamfile
@@ -3,8 +3,8 @@
# Optimization and Debug flags
-PREF_CCFLAGS += $(CCOPTFLAG) ; # Turn optimisation on
-#PREF_CCFLAGS += $(CCDEBUGFLAG) ; # Debugging flags
+#PREF_CCFLAGS += $(CCOPTFLAG) ; # Turn optimisation on
+PREF_CCFLAGS += $(CCDEBUGFLAG) ; # Debugging flags
#PREF_CCFLAGS += $(CCPROFFLAG) ; # Profile flags
#PREF_LINKFLAGS += $(LINKPROFFLAG) ; # Profile flags
#PREF_CCFLAGS += $(CCHEAPDEBUG) ; # Heap Debugging flags
@@ -27,9 +27,9 @@ HDRS = ../h ../numlib ../plot $(TIFFINC) ;
LINKLIBS = librspl ../numlib/libnum ../plot/libplot ../plot/libvrml ../icc/libicc $(TIFFLIB) $(JPEGLIB) ;
# Test programs
-MainsFromSources revbench.c c1.c c1df.c t2d.c t2ddf.c t3d.c t3ddf.c tnd.c trnd.c ;
+MainsFromSources revbench.c c1.c cw1.c cw3.c c1df.c t2d.c t2ddf.c t3d.c t3ddf.c tnd.c trnd.c ;
-BUILD_TESTS = true ;
+BUILD_TESTS = false ;
if $(BUILD_TESTS) {
diff --git a/rspl/Makefile.am b/rspl/Makefile.am
deleted file mode 100644
index cfa0ec2..0000000
--- a/rspl/Makefile.am
+++ /dev/null
@@ -1,15 +0,0 @@
-include $(top_srcdir)/Makefile.shared
-
-privatelib_LTLIBRARIES = librspl.la
-privatelibdir = $(pkglibdir)
-
-librspl_la_SOURCES = rspl.h rspl_imp.h mlbs.h rspl.c scat.c rev.c \
- rev.h gam.c spline.c opt.c
-librspl_la_LIBADD = ../numlib/libargyllnum.la ../plot/libvrml.la
-
-LDADD = ./librspl.la ../numlib/libargyllnum.la ../plot/libplot.la \
- ../plot/libvrml.la $(X_LIBS) $(TIFF_LIBS) $(ICC_LIBS)
-
-check_PROGRAMS = revbench c1 c1df t2d t2ddf t3d t3ddf tnd trnd
-
-EXTRA_DIST = License.txt Readme.txt
diff --git a/rspl/afiles b/rspl/afiles
index af20c44..7d9c29f 100644
--- a/rspl/afiles
+++ b/rspl/afiles
@@ -16,6 +16,8 @@ gam.h
gam.c
revbench.c
c1.c
+cw1.c
+cw3.c
c1df.c
spline.c
t2d.c
diff --git a/rspl/c1.c b/rspl/c1.c
index 8172cef..2717920 100644
--- a/rspl/c1.c
+++ b/rspl/c1.c
@@ -19,7 +19,7 @@
#undef DIAG
#undef DIAG2
#undef GLOB_CHECK
-#undef RES2 /* Do multiple test at various resolutions */
+#define RES2 /* Do multiple test at various resolutions */
#define AVGDEV 0.0 /* Average deviation of function data */
#include <stdio.h>
@@ -285,7 +285,6 @@ int main(int argc, char *argv[]) {
rss->fit_rspl(rss,
0 | (twopass ? RSPL_2PASSSMTH : 0)
| (extra ? RSPL_EXTRAFIT2 : 0) ,
- 0,
test_points, /* Test points */
pnts, /* Number of test points */
low, high, gres, /* Low, high, resolution of grid */
diff --git a/rspl/cw1.c b/rspl/cw1.c
new file mode 100644
index 0000000..75a4f55
--- /dev/null
+++ b/rspl/cw1.c
@@ -0,0 +1,308 @@
+
+/************************************************/
+/* Investigate various curve approximations */
+/************************************************/
+
+/* Discrete regularized spline versions */
+/* Test variable grid spacing in 1D
+
+/* Author: Graeme Gill
+ * Date: 4/10/95
+ * Date: 5/4/96
+ *
+ * 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 "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 10
+#define GRES1 100
+
+/* 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.5);
+
+ if (sgn)
+ out = -out;
+
+ return out;
+}
+
+/* Grid spacing power */
+#define GRIDPOW 1.5
+
+/* Scattered data sample spacing power */
+#define SAMPPOW 1.3
+
+
+co test_points[PNTS1];
+
+double ipos[GRES1];
+double *iposes[1] = { ipos };
+
+double powlike(double vv, double pp);
+
+void usage(void) {
+ fprintf(stderr,"Test 1D 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, 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;
+ high[0] = 1.0;
+ avgdev[0] = 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;
+ }
+
+ /* Only one trial */
+ for (n = 0; n < TRIALS; n++) {
+ double lrand = 0.0; /* Amount of level randomness */
+ int pnts;
+ int fres;
+
+ pnts = PNTS1;
+ fres = GRES1;
+ gres[0] = fres;
+
+ /* Create the object */
+ rss = new_rspl(RSPL_NOFLAGS,
+ 1, /* di */
+ 1); /* fdi */
+
+ /* For each grid spacing test */
+ for (j = 0; j < 2; j++) {
+
+ for (i = 0; i < pnts; i++) {
+ double val, out;
+
+ val = i/(pnts - 1.0);
+
+ if (j == 0) { /* Linear */
+// out = func(val + d_rand(-0.001, 0.001));
+ out = func(val);
+ test_points[i].p[0] = val;
+ test_points[i].v[0] = 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[0] = val;
+ test_points[i].v[0] = out;
+ }
+ }
+
+ /* Set grid position info */
+ 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;
+ }
+
+ rss->fit_rspl(rss,
+ 0,
+ test_points, /* Test points */
+ pnts, /* 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 */
+ 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;
+ 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;
+ 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 */
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 */
diff --git a/rspl/rev.c b/rspl/rev.c
index 237dbe3..1dbf6cc 100644
--- a/rspl/rev.c
+++ b/rspl/rev.c
@@ -22,11 +22,24 @@
Should fix the clipping case so that a direction weighting
funtion can be applied. This should be used just like
the perceptual case to increase L* constance for dark
- colors. Will this stuff up the geometric consistency though ?
- [ See fill_nncell(), fill_nncell() and users of calc_fwd_nn_cell_list(),
- ie. nnearest_clip_solve(), clipn_setsort() etc. ]
+ colors. This would entail large scale changes though,
+ since a lot of code assumes minimal euclidean distance
+ goal, from the cell selection structure [ See fill_nncell(),
+ fill_nncell() and users of calc_fwd_nn_cell_list() ] and
+ the within cell computation [ ie. See nnearest_clip_solve(),
+ clipn_setsort() etc. ]
+ XYZ PCS couldn't work with a simple weighting - it would have
+ to be a position dependent weighting.
The SVD least squares computation case makes this hard to change ?
Would have to feed in a weighting function, or can it be general ?
+ -
+ Can this be solved some other way, ie. by using gamut
+ mapping type look up ? Problem is precision.
+ -
+ Vector clip could be used (if intent can be turned
+ into computable vector clip direction), but it is slow,
+ because it search all cells from source until it
+ hits surface.
Allow function callback to set auxiliary values for
flag RSPL_AUXLOCUS.
@@ -1437,7 +1450,8 @@ unsigned int tcount /* grid touch count for this operation */
nunlk++;
}
fprintf(stdout,"Diagnostic: rev.sz = %lu, rev.max_sz = %lu, numlocked = %d, nunlk = %d\n",
- rc->s->rev.sz, rc->s->rev.max_sz, rc->nunlocked,nunlk);
+ (unsigned long)rc->s->rev.sz, (unsigned long)rc->s->rev.max_sz,
+ rc->nunlocked, nunlk);
error("Not enough memory to process in chunks");
}
break; /* cache has run out of room, so abandon, and do it next time */
@@ -4971,7 +4985,7 @@ rspl *s /* Pointer to rspl grid */
g_no_rev_cache_instances > 1 ? "are" : "is",
g_no_rev_cache_instances,
g_no_rev_cache_instances > 1 ? "s" : "",
- ram_portion/1000000);
+ (unsigned long)ram_portion/1000000);
}
}
@@ -6165,7 +6179,7 @@ if (prop != NULL) {
g_no_rev_cache_instances > 1 ? "are" : "is",
g_no_rev_cache_instances,
g_no_rev_cache_instances > 1 ? "s" : "",
- ram_portion/1000000);
+ (unsigned long)ram_portion/1000000);
}
s->rev.rev_valid = 1;
@@ -6233,7 +6247,7 @@ rspl *s /* Pointer to rspl grid */
g_no_rev_cache_instances > 1 ? "are" : "is",
g_no_rev_cache_instances,
g_no_rev_cache_instances > 1 ? "s" : "",
- ram_portion/1000000);
+ (unsigned long)ram_portion/1000000);
}
}
s->rev.rev_valid = 0;
@@ -6497,7 +6511,7 @@ rspl *s
if (g_avail_ram > safe_max_vmem) {
g_avail_ram = safe_max_vmem;
if (s->verbose && repsr == 0)
- fprintf(stdout,"%cTrimmed maximum cache RAM to %lu Mbytes to allow for VM limit\n",cr_char,g_avail_ram/1000000);
+ fprintf(stdout,"%cTrimmed maximum cache RAM to %lu Mbytes to allow for VM limit\n",cr_char,(unsigned long)g_avail_ram/1000000);
}
}
@@ -6516,7 +6530,7 @@ rspl *s
}
if (max_vmem != 0 && g_avail_ram > max_vmem && repsr == 0) {
g_avail_ram = (size_t)(0.95 * max_vmem);
- fprintf(stdout,"%cARGYLL_REV_CACHE_MULT * RAM trimmed to %lu Mbytes to allow for VM limit\n",cr_char,g_avail_ram/1000000);
+ fprintf(stdout,"%cARGYLL_REV_CACHE_MULT * RAM trimmed to %lu Mbytes to allow for VM limit\n",cr_char,(unsigned long)g_avail_ram/1000000);
}
}
@@ -6525,7 +6539,7 @@ rspl *s
DBG(("reverse cache max memory = %d Mbytes\n",s->rev.max_sz/1000000));
if (s->verbose && repsr == 0) {
- fprintf(stdout, "%cRev cache RAM = %lu Mbytes\n",cr_char,g_avail_ram/1000000);
+ fprintf(stdout, "%cRev cache RAM = %lu Mbytes\n",cr_char,(unsigned long)g_avail_ram/1000000);
repsr = 1;
}
diff --git a/rspl/rspl.c b/rspl/rspl.c
index 13a1776..23ef6b4 100644
--- a/rspl/rspl.c
+++ b/rspl/rspl.c
@@ -866,7 +866,7 @@ rspl *s
/* ============================================ */
/* Initialize the grid from a provided function. By default the grid */
-/* values are set to exactly the value returned fy func(), unless the */
+/* values are set to exactly the value returned by func(), unless the */
/* RSPL_SET_APXLS flag is set, in which case an attempt is made to have */
/* the grid points represent a least squares aproximation to the underlying */
/* surface, by using extra samples in the middle of grid cells. */
diff --git a/rspl/scat.c b/rspl/scat.c
index b4ed978..ec7e99c 100644
--- a/rspl/scat.c
+++ b/rspl/scat.c
@@ -128,10 +128,6 @@
/* algorithm parameters [Release defaults] */
#undef POINTWEIGHT /* [Undef] Increas smoothness weighting proportional to number of points */
#define INCURVEADJ /* [Defined] Adjust smoothness criteria for input curve grid spacing */
-#define EXTRA_SURFACE_SMOOTHING /* [Defined] Stiffen surface points to comp. for single ended. */
- /* The following are available, but the smoothing table is */
- /* not setup for them, and they are not sufficiently different */
- /* from the default smoothing to be useful. */
#define ENABLE_2PASSSMTH /* [Define] Enable 2 pass smooth using Gaussian filter */
#define ENABLE_EXTRAFIT /* [Undef] Enable the extra fit option. Good to combat high smoothness. */
@@ -142,12 +138,12 @@
/* Experimental set: */
-#pragma message("!!!!!!!!! Experimental config set !!!!!!!!!")
+#pragma message("!!!!!!!!! Experimental hi-accuracy config set !!!!!!!!!")
#define TOL 1e-12 /* Tollerance of result - usually 1e-5 is best. */
#define TOL_IMP 1.0 /* Minimum error improvement to continue - reduces accuracy (1.0 == off) */
#undef GRADUATED_TOL /* Speedup attemp - use reduced tollerance for prior grids. */
-#define GRATIO 2.0 /* Multi-grid resolution ratio */
+#define GRATIO 1.5 /* Multi-grid resolution ratio */
#undef OVERRLX /* Use over relaxation factor when progress slows (worse accuracy ?) */
#define JITTERS 0 /* Number of 1D conjugate solve itters */
#define CONJ_TOL 1.0 /* Extra tolereance on 1D conjugate solution times TOL. */
@@ -160,8 +156,8 @@
/* Release set: */
-#define TOL 1e-6 /* [1e-6] Tollerance of result - usually 1e-5 is best. */
-#define TOL_IMP 0.998 /* [0.998] Minimum error improvement to continue - reduces accuracy (1.0 == off) */
+#define TOL 1e-7 /* [1e-6] Tollerance of result - usually 1e-5 is best. */
+#define TOL_IMP 0.999 /* [0.998] Minimum error improvement to continue - reduces accuracy (1.0 == off) */
#undef GRADUATED_TOL /* [Undef] Speedup attemp - use reduced tollerance for prior grids. */
#define GRATIO 2.0 /* [2.0] Multi-grid resolution ratio */
#undef OVERRLX /* [Undef] Use over relaxation when progress slows (worse accuracy ?) */
@@ -262,7 +258,7 @@ static void free_cj_arrays(cj_arrays *ta);
static int add_rspl_imp(rspl *s, int flags, void *d, int dtp, int dno);
static mgtmp *new_mgtmp(rspl *s, int gres[MXDI], int f);
static void free_mgtmp(mgtmp *m);
-static void setup_solve(mgtmp *m, int final);
+static void setup_solve(mgtmp *m);
static void solve_gres(mgtmp *m, cj_arrays *ta, double tol, int final);
static void init_soln(mgtmp *m1, mgtmp *m2);
static void comp_ccv(mgtmp *m);
@@ -656,8 +652,7 @@ add_rspl_imp(
if (s->tpsm && s->tpsm2 != 0) { /* 2nd pass of 2 pass smoothing */
init_ccv(m); /* Downsample m->ccv from s->g.ccv */
}
-// setup_solve(m, nn == (s->niters-1));
- setup_solve(m, 1);
+ setup_solve(m);
if (nn == 0) { /* Make sure we have an initial x[] */
for (i = 0; i < m->g.no; i++)
@@ -921,9 +916,20 @@ free_data(rspl *s) {
/* to them, would have an average sample deviation of 0.05. */
/* For normally distributed errors, the average deviation is */
/* aproximately 0.564 times the standard deviation. (0.564 * sqrt(variance)) */
-/* This table is appropriate for the default rspl algorithm + slight EXTRA_SURFACE_SMOOTHING, */
+/* This table is appropriate for the default rspl algorithm, */
/* and is NOT setup for RSPL_2PASSSMTH or RSPL_EXTRAFIT2 !! */
/* SMOOTH */
+
+/* There are still issues with all this - the level of smoothing actually */
+/* depends on the degree of fit of the underlying model - ie. how close */
+/* to straight the mapping is. To get actual noise reduction under these */
+/* conditions is harder than when there is some curvature to "tension" things. */
+/* This is evident is Lab vs. XYZ display profiles, and there is code */
+/* in xlut.c that tries to adjust to this. */
+
+/* !!! Possible answer - should be using third order differences */
+/* for controlling smoothness, not second order (curvature) ?? */
+
// ~~99
static double opt_smooth(
rspl *s,
@@ -963,7 +969,7 @@ static double opt_smooth(
/* New for V1.10, from smtmpp using sRGB, EpsonR1800, Hitachi2112, */
- /* Fogra39L, Canon1180, Epson10K, with low EXTRA_SURFACE_SMOOTHING. */
+ /* Fogra39L, Canon1180, Epson10K (did use EXTRA_SURFACE_SMOOTHING). */
/* Main lookup table, by [di][ncix][adix]: */
/* Values are log of smoothness value. */
@@ -1349,8 +1355,7 @@ static void free_mgtmp(mgtmp *m) {
*/
static void setup_solve(
- mgtmp *m, /* initialized grid temp structure */
- int final /* nz if final resolution (activate EXTRA_SURFACE_SMOOTHING) */
+ mgtmp *m /* initialized grid temp structure */
) {
rspl *s = m->s;
int di = s->di;
@@ -1470,7 +1475,6 @@ static void setup_solve(
b[i] = 0.0;
}
-#ifdef NEVER
/* Production version, without extra edge weight */
/* Accumulate curvature dependent factors to the triangular A matrix. */
@@ -1499,20 +1503,46 @@ static void setup_solve(
ki-1 * w[i-1] * w[i-1] * u[i-1]
ki * -(w[i-1] + w[i]) * -(w[i-1] + w[i]) * u[i])
ki+1 * w[i] * w[i] * u[i+1]
+
+ The weighting is adjusted to compensate for increased
+ grid density by the dw weight, and for the different
+ scaling either side of a grid point with the w0 and w1 scale.
+ Both are normalised so that the overall smoothness scaling
+ factor applies.
*/
- { /* Setting this up from scratch */
+ { /* Setting this up from scratch */
+ double dwtw[MXDIDO]; /* Density weight normalizer */
+ int wtwc[MXDIDO];
ECOUNT(gc, MXDIDO, di, 0, gres, 0);
- EC_INIT(gc);
for (oawt = 0.0, i = 1; i < 21; i++)
oawt += wvals[i];
oawt *= wvals[0];
+#define EN_DW /* Enable grid density weighting */
+#define EN_CW /* Enable grid scale weighting */
+
+ /* Compute the ipos[] weight normalisation factors */
+ for (e = 0; e < di; e++) {
+ if (m->g.ipos[e] == NULL)
+ continue;
+ dwtw[e] = 0.0;
+ for (i = 1; i < gres[e]; i++) {
+ double w;
+ w = fabs(m->g.ipos[e][i] - m->g.ipos[e][i-1]);
+//printf("[%d][%d] w = %f\n",e,i,w);
+ dwtw[e] += w;
+ }
+ dwtw[e] /= (gres[e] - 1.0); /* Average weights */
+//printf("dwtw[%d] = %f\n",e,dwtw[e]);
+ }
+
+ EC_INIT(gc);
for (i = 0; i < gno; i++) {
for (e = 0; e < di; e++) { /* For each curvature direction */
- double w0, w1, tt;
+ double dw, w0, w1, tt;
double cw = 2.0 * m->sf.cw[e]; /* Overall curvature weight */
cw *= s->d.vw[f]; /* Scale curvature weight for data range */
@@ -1521,14 +1551,24 @@ static void setup_solve(
if ((gc[e]-2) >= 0) {
/* double kw = cw * gp[UO_C(e,1)].k; */ /* Cell bellow k value */
double kw = cw;
- w0 = w1 = 1.0;
+ dw = w1 = 1.0;
if (m->g.ipos[e] != NULL) {
+
w0 = fabs(m->g.ipos[e][gc[e]-1] - m->g.ipos[e][gc[e]-2]);
w1 = fabs(m->g.ipos[e][gc[e]-0] - m->g.ipos[e][gc[e]-1]);
- tt = sqrt(w0 * w1); /* Normalise overall width weighting effect */
+//printf("raw [%d][%d] w0 = %f, w1 = %f\n",gc[e],i,w0,w1);
+#ifdef EN_DW
+ dw = 0.5 * (w0 + w1); /* in this direction around -1 */
+ dw = dw/dwtw[e]; /* normalise */
+#endif
+ tt = sqrt(w0 * w1);
w1 = tt/w1;
+#ifndef EN_CW
+ w1 = 1.0;
+#endif
+//printf("[%d][%d] dw = %f, w1 = %f\n",gc[e],i,dw,w1);
}
- A[i][ixcol[0]] += w1 * w1 * kw;
+ A[i][ixcol[0]] += dw * w1 * w1 * kw;
if (ccv != NULL)
b[i] += kw * (w1) * ccv[i - gci[e]][e]; /* Curvature compensation value */
}
@@ -1537,16 +1577,25 @@ static void setup_solve(
if ((gc[e]-1) >= 0 && (gc[e]+1) < gres[e]) {
/* double kw = cw * gp->k; */ /* This cells k value */
double kw = cw;
- w0 = w1 = 1.0;
+ dw = w0 = w1 = 1.0;
if (m->g.ipos[e] != NULL) {
w0 = fabs(m->g.ipos[e][gc[e]-0] - m->g.ipos[e][gc[e]-1]);
w1 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]-0]);
+//printf("raw [%d][%d] w0 = %f, w1 = %f\n",gc[e],i,w0,w1);
+#ifdef EN_DW
+ dw = 0.5 * (w0 + w1); /* in this direction around 0 */
+ dw = dw/dwtw[e]; /* normalise */
+#endif
tt = sqrt(w0 * w1);
w0 = tt/w0;
w1 = tt/w1;
+#ifndef EN_CW
+ w0 = w1 = 1.0;
+#endif
+//printf("[%d][%d] dw = %f, w0 = %f, w1 = %f\n",gc[e],i,dw,w0,w1);
}
- A[i][ixcol[0]] += -(w0 + w1) * -(w0 + w1) * kw;
- A[i][ixcol[gci[e]]] += -(w0 + w1) * w1 * kw * oawt;
+ A[i][ixcol[0]] += dw * -(w0 + w1) * -(w0 + w1) * kw;
+ A[i][ixcol[gci[e]]] += dw * -(w0 + w1) * w1 * kw * oawt;
if (ccv != NULL)
b[i] += kw * -(w0 + w1) * ccv[i][e]; /* Curvature compensation value */
}
@@ -1555,173 +1604,33 @@ static void setup_solve(
if ((gc[e]+2) < gres[e]) {
/* double kw = cw * gp[UO_C(e,2)].k; */ /* Cell above k value */
double kw = cw;
- w0 = w1 = 1.0;
+ dw = w0 = w1 = 1.0;
if (m->g.ipos[e] != NULL) {
w0 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]+0]);
w1 = fabs(m->g.ipos[e][gc[e]+2] - m->g.ipos[e][gc[e]+1]);
- tt = sqrt(w0 * w1);
- w0 = tt/w0;
- w1 = tt/w1;
- }
- A[i][ixcol[0]] += w0 * w0 * kw;
- A[i][ixcol[gci[e]]] += w0 * -(w0 + w1) * kw;
- A[i][ixcol[2 * gci[e]]] += w0 * w1 * kw;
- if (ccv != NULL)
- b[i] += kw * -(w0 + w1) * ccv[i][e]; /* Curvature compensation value */
- }
- }
- EC_INC(gc);
- }
- }
-#endif /* NEVER */
-
-#ifdef ALWAYS
- /* Production version that allows for extra weight on grid edges */
-
- /* Accumulate curvature dependent factors to the triangular A matrix. */
- /* Because it's triangular, we compute and add in all the weighting */
- /* factors at and to the right of each cell. */
-
- /* The ipos[] factor is to allow for the possibility that the */
- /* grid spacing may be non-uniform in the colorspace where the */
- /* function being modelled is smooth. Our curvature computation */
- /* needs to make allowsance for this fact in computing the */
- /* node value differences that equate to zero curvature. */
- /*
- The old curvature fixed grid spacing equation was:
- ki * (u[i-1] - 2 * u[i] + u[i+1])^2
- with derivatives wrt each node:
- ki-1 * 1 * 2 * u[i-1]
- ki * -2 * 2 * u[i]
- ki+1 * 1 * 2 * u[i+1]
-
- Allowing for scaling of each grid difference by w[i-1] and w[i],
- where w[i-1] corresponds to the width of cell i-1 to i,
- and w[i] corresponds to the width of cell i to i+1:
- ki * (w[i-1] * (u[i-1] - u[i]) + w[i] * (u[i+1] - u[i[))^2
- = ki * (w[i-1] * u[i-1] - (w[i-1] + w[i]) * u[i]) + w[i] * u[i+1])^2
- with derivatives wrt each node:
- ki-1 * w[i-1] * w[i-1] * u[i-1]
- ki * -(w[i-1] + w[i]) * -(w[i-1] + w[i]) * u[i])
- ki+1 * w[i] * w[i] * u[i+1]
- */
- { /* Setting this up from scratch */
- ECOUNT(gc, MXDIDO, di, 0, gres, 0);
-#ifdef EXTRA_SURFACE_SMOOTHING
-// double k0w = 4.0, k1w = 1.3333; /* Extra stiffness */
-// double k0w = 3.0, k1w = 1.26; /* Some extra stiffness */
- double k0w = 2.0, k1w = 1.15; /* A little extra stiffness */
-#else
- double k0w = 1.0, k1w = 1.0; /* No extra weights */
+//printf("raw [%d][%d] w0 = %f, w1 = %f\n",gc[e],i,w0,w1);
+#ifdef EN_DW
+ dw = 0.5 * (w0 + w1); /* in this direction around +1 */
+ dw = dw/dwtw[e]; /* normalise */
#endif
-
- EC_INIT(gc);
- for (oawt = 0.0, i = 1; i < 21; i++)
- oawt += wvals[i];
- oawt *= wvals[0];
-
- if (final == 0)
- k0w = k1w = 1.0; /* Activate extra edge smoothing on final grid ? */
-
- for (i = 0; i < gno; i++) {
- int k;
-
- /* We're creating the equation cooeficients for solving the */
- /* partial derivative equation w.r.t. node point i. */
- /* Due to symetry in the smoothness interactions, only */
- /* the triangle cooeficients of neighbour nodes is needed. */
- for (e = 0; e < di; e++) { /* For each curvature direction */
- double kw, w0, w1, tt;
- double cw = 2.0 * m->sf.cw[e]; /* Overall curvature weight */
- double xx = 1.0; /* Extra edge weighing */
- cw *= s->d.vw[f]; /* Scale curvature weight for data range */
-
- /* weight factor for outer or 2nd outer in other dimensions */
- for (k = 0; k < di; k++) {
- if (k == e)
- continue;
- if (gc[k] == 0 || gc[k] == (gres[k]-1))
- xx *= k0w;
- else if (gc[k] == 1 || gc[k] == (gres[k]-2))
- xx *= k1w;
- }
-
- /* If at least two above lower edge in this dimension */
- /* Add influence on Curvature of cell below */
- if ((gc[e]-2) >= 0) {
- /* double kw = cw * gp[-gc[e]].k; */ /* Cell bellow k value */
- kw = cw * xx;
- w0 = w1 = 1.0;
- if (m->g.ipos[e] != NULL) {
- w0 = fabs(m->g.ipos[e][gc[e]-1] - m->g.ipos[e][gc[e]-2]);
- w1 = fabs(m->g.ipos[e][gc[e]-0] - m->g.ipos[e][gc[e]-1]);
- tt = sqrt(w0 * w1); /* Normalise overall width weighting effect */
- w1 = tt/w1;
- }
- if ((gc[e]-2) == 0 || (gc[e]+0) == (gres[e]-1))
- kw *= k0w;
- else if ((gc[e]-2) == 1 || (gc[e]-0) == (gres[e]-2))
- kw *= k1w;
- A[i][ixcol[0]] += kw * (w1) * w1;
- if (ccv != NULL) {
-//printf("~1 tweak b[%d] by %e\n",i,kw * (w1) * ccv[i - gci[e]][e]);
- b[i] += kw * (w1) * ccv[i - gci[e]][e]; /* Curvature compensation value */
- }
- }
- /* If not one from upper or lower edge in this dimension */
- /* Add influence on Curvature of this cell */
- if ((gc[e]-1) >= 0 && (gc[e]+1) < gres[e]) {
- /* double kw = cw * gp->k; */ /* This cells k value */
- kw = cw * xx;
- w0 = w1 = 1.0;
- if (m->g.ipos[e] != NULL) {
- w0 = fabs(m->g.ipos[e][gc[e]-0] - m->g.ipos[e][gc[e]-1]);
- w1 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]-0]);
tt = sqrt(w0 * w1);
w0 = tt/w0;
w1 = tt/w1;
+#ifndef EN_CW
+ w0 = w1 = 1.0;
+#endif
+//printf("[%d][%d] dw = %f, w0 = %f, w1 = %f\n",gc[e],i,dw,w0,w1);
}
- if ((gc[e]-1) == 0 || (gc[e]+1) == (gres[e]-1))
- kw *= k0w;
- else if ((gc[e]-1) == 1 || (gc[e]+1) == (gres[e]-2))
- kw *= k1w;
- A[i][ixcol[0]] += kw * -(w0 + w1) * -(w0 + w1);
- A[i][ixcol[gci[e]]] += kw * -(w0 + w1) * w1 * oawt;
- if (ccv != NULL) {
-//printf("~1 tweak b[%d] by %e\n",i, kw * -(w0 + w1) * ccv[i][e]);
+ A[i][ixcol[0]] += dw * w0 * w0 * kw;
+ A[i][ixcol[gci[e]]] += dw * w0 * -(w0 + w1) * kw;
+ A[i][ixcol[2 * gci[e]]] += dw * w0 * w1 * kw;
+ if (ccv != NULL)
b[i] += kw * -(w0 + w1) * ccv[i][e]; /* Curvature compensation value */
- }
- }
- /* If at least two below the upper edge in this dimension */
- /* Add influence on Curvature of cell above */
- if ((gc[e]+2) < gres[e]) {
- /* double kw = cw * gp[gc[e]].k; */ /* Cell above k value */
- kw = cw * xx;
- w0 = w1 = 1.0;
- if (m->g.ipos[e] != NULL) {
- w0 = fabs(m->g.ipos[e][gc[e]+1] - m->g.ipos[e][gc[e]+0]);
- w1 = fabs(m->g.ipos[e][gc[e]+2] - m->g.ipos[e][gc[e]+1]);
- tt = sqrt(w0 * w1);
- w0 = tt/w0;
- w1 = tt/w1;
- }
- if ((gc[e]+0) == 0 || (gc[e]+2) == (gres[e]-1))
- kw *= k0w;
- else if ((gc[e]+0) == 1 || (gc[e]+2) == (gres[e]-2))
- kw *= k1w;
- A[i][ixcol[0]] += kw * (w0) * w0;
- A[i][ixcol[gci[e]]] += kw * (w0) * -(w0 + w1);
- A[i][ixcol[2 * gci[e]]] += kw * (w0) * w1;
- if (ccv != NULL) {
-//printf("~1 tweak b[%d] by %e\n",i, kw * (w0) * ccv[i + gci[e]][e]);
- b[i] += kw * (w0) * ccv[i + gci[e]][e]; /* Curvature compensation value */
- }
}
}
EC_INC(gc);
}
}
-#endif /* ALWAYS */
#ifdef DEBUG
printf("After adding curvature equations:\n");
@@ -1736,7 +1645,6 @@ static void setup_solve(
nbsum = 0.0; /* Zero sum of b[] squared */
-#ifdef ALWAYS
/* Accumulate weak default function factors. These are effectively a */
/* weak "data point" exactly at each grid point. */
/* (Note we're not currently doing this in a cache friendly order, */
@@ -1775,9 +1683,7 @@ static void setup_solve(
#endif /* DEBUG */
}
-#endif /* ALWAYS */
-#ifdef ALWAYS
/* Accumulate data point dependent factors */
for (n = 0; n < dno; n++) { /* Go through all the data points */
int j,k;
@@ -1815,8 +1721,6 @@ static void setup_solve(
nbsum = 1e-4;
m->q.normb = nbsum;
-#endif /* ALWAYS */
-
#ifdef DEBUG
printf("After adding data point equations:\n");
for (i = 0; i < gno; i++) {