summaryrefslogtreecommitdiff
path: root/rspl
diff options
context:
space:
mode:
Diffstat (limited to 'rspl')
-rw-r--r--rspl/Jamfile2
-rwxr-xr-x[-rw-r--r--]rspl/c1.c0
-rwxr-xr-x[-rw-r--r--]rspl/c1df.c0
-rwxr-xr-x[-rw-r--r--]rspl/cw1.c0
-rwxr-xr-x[-rw-r--r--]rspl/cw3.c0
-rwxr-xr-x[-rw-r--r--]rspl/gam.c13
-rwxr-xr-x[-rw-r--r--]rspl/gam.h0
-rwxr-xr-x[-rw-r--r--]rspl/mlbs.c0
-rwxr-xr-x[-rw-r--r--]rspl/mlbs.h0
-rwxr-xr-x[-rw-r--r--]rspl/opt.c0
-rwxr-xr-x[-rw-r--r--]rspl/rev.c963
-rwxr-xr-x[-rw-r--r--]rspl/rev.h17
-rwxr-xr-x[-rw-r--r--]rspl/revbench.c0
-rwxr-xr-x[-rw-r--r--]rspl/rspl.c0
-rwxr-xr-x[-rw-r--r--]rspl/rspl.h35
-rwxr-xr-x[-rw-r--r--]rspl/rspl1.c0
-rwxr-xr-x[-rw-r--r--]rspl/rspl1.h0
-rwxr-xr-x[-rw-r--r--]rspl/rspl_imp.h0
-rwxr-xr-x[-rw-r--r--]rspl/scat.c127
-rwxr-xr-x[-rw-r--r--]rspl/scat2.c0
-rwxr-xr-x[-rw-r--r--]rspl/sm1.c0
-rwxr-xr-x[-rw-r--r--]rspl/sm2.c0
-rwxr-xr-x[-rw-r--r--]rspl/sm3.c0
-rwxr-xr-x[-rw-r--r--]rspl/smtmpp.c0
-rwxr-xr-x[-rw-r--r--]rspl/smtnd.c2
-rwxr-xr-x[-rw-r--r--]rspl/spline.c0
-rwxr-xr-x[-rw-r--r--]rspl/stest.c0
-rwxr-xr-x[-rw-r--r--]rspl/t2d.c0
-rwxr-xr-x[-rw-r--r--]rspl/t2ddf.c0
-rwxr-xr-x[-rw-r--r--]rspl/t3d.c0
-rwxr-xr-x[-rw-r--r--]rspl/t3ddf.c0
-rwxr-xr-x[-rw-r--r--]rspl/tnd.c0
-rwxr-xr-x[-rw-r--r--]rspl/trnd.c0
33 files changed, 917 insertions, 242 deletions
diff --git a/rspl/Jamfile b/rspl/Jamfile
index 0d93acc..6f80f9f 100644
--- a/rspl/Jamfile
+++ b/rspl/Jamfile
@@ -27,6 +27,8 @@ HDRS = ../h ../numlib ../plot $(TIFFINC) ;
LINKLIBS = librspl ../plot/libplot ../numlib/libnum ../numlib/libui ../plot/libvrml ../icc/libicc $(TIFFLIB) $(JPEGLIB) ;
# Test programs
+LINKFLAGS += $(GUILINKFLAGS) ;
+
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 = false ;
diff --git a/rspl/c1.c b/rspl/c1.c
index 31849ee..31849ee 100644..100755
--- a/rspl/c1.c
+++ b/rspl/c1.c
diff --git a/rspl/c1df.c b/rspl/c1df.c
index 4020d7a..4020d7a 100644..100755
--- a/rspl/c1df.c
+++ b/rspl/c1df.c
diff --git a/rspl/cw1.c b/rspl/cw1.c
index 2a026c8..2a026c8 100644..100755
--- a/rspl/cw1.c
+++ b/rspl/cw1.c
diff --git a/rspl/cw3.c b/rspl/cw3.c
index e26cb7a..e26cb7a 100644..100755
--- a/rspl/cw3.c
+++ b/rspl/cw3.c
diff --git a/rspl/gam.c b/rspl/gam.c
index 0765e05..2a9a63b 100644..100755
--- a/rspl/gam.c
+++ b/rspl/gam.c
@@ -21,9 +21,14 @@
/*
This probably needs re-writing, since (I think) it doesn't take gamut
- self intersection into account. Outline would be:
+ self intersection into account. The rspl/rec.c code works more
+ reliably, but doesn't cope with concavities that cross the line
+ to the focal point. So can this "surface following" code be
+ adapted to solve this type of problem ?
- Have a spactial cache structure that contains list of potentialy
+ Outline would be:
+
+ Have a spatial cache structure that contains list of potentialy
intersecting triangles for any point. Create this incrementally.
Each triangle can be replaced by a list of fragment triangles
sharing the same plane.
@@ -61,6 +66,9 @@
Each intersection adds new nodes, and splits a triangle into
about 8 smaller triangles. Trick is to avoid slivers and
numerical issues with triangles.
+
+ Alternative would be to not split triangles, but just track
+ the line of intersection ?
*/
@@ -72,6 +80,7 @@
Add ink limit support. This be done by breaking
a cell into a fixed geometry of smaller simplexes
by dividing the cell into two on each axis.
+ (rspl/rev.c has code for this ??)
Need to then add scan that detects areas to prune,
that then ties in with rev code to mark such
diff --git a/rspl/gam.h b/rspl/gam.h
index 0892e4d..0892e4d 100644..100755
--- a/rspl/gam.h
+++ b/rspl/gam.h
diff --git a/rspl/mlbs.c b/rspl/mlbs.c
index bbe3865..bbe3865 100644..100755
--- a/rspl/mlbs.c
+++ b/rspl/mlbs.c
diff --git a/rspl/mlbs.h b/rspl/mlbs.h
index 4678cfd..4678cfd 100644..100755
--- a/rspl/mlbs.h
+++ b/rspl/mlbs.h
diff --git a/rspl/opt.c b/rspl/opt.c
index d5a70ce..d5a70ce 100644..100755
--- a/rspl/opt.c
+++ b/rspl/opt.c
diff --git a/rspl/rev.c b/rspl/rev.c
index f2f717b..017943b 100644..100755
--- a/rspl/rev.c
+++ b/rspl/rev.c
@@ -72,7 +72,6 @@
a particular aux target that lies within that segment.
(1150 near black, k ~= 0.4).
-
*/
#include <stdio.h>
@@ -109,9 +108,13 @@
//#undef DMALLOC_GLOBALS
#define DOSORT /* [def] Cell sort for better speed */
+#undef EN_UNTWIST /* [und] Force attempt to try and untwist gamut surface */
+ /* - Seems to improve some, make some worse ?? (i.e Bonet) */
+ /* By default this is controlled using ARGYLL_UNTWIST_GAMUT_SURFACE */
+ /* environment variable. */
#undef REVTABLESTATS /* [und] Reverse table stats */
-#undef REVVRML /* [und] Reverse table plots */
+#undef REVVRML /* [und] Reverse table 3D plots */
#undef DEBUG1 /* [und] Higher level code */
#undef DEBUG2 /* [und] Lower level code */
@@ -172,16 +175,16 @@ int thissz, lastsz = -1;
#endif
#ifdef REVTABLESTATS
-#pragma message("!!!!!!!!! REVTABLESTATS set in rspl/rev.c !!!!!!!!!!!")
+# pragma message("!!!!!!!!! REVTABLESTATS set in rspl/rev.c !!!!!!!!!!!")
#endif
#ifdef REVVRML
-#pragma message("!!!!!!!!! REVVRML set in rspl/rev.c !!!!!!!!!!!")
-#include "vrml.h"
+# pragma message("!!!!!!!!! REVVRML set in rspl/rev.c !!!!!!!!!!!")
+# include "vrml.h"
#endif
#ifdef CHECK_NNLU
-#pragma message("!!!!!!!!! CHECK_NNLU set in rspl/rspl.h !!!!!!!!!!!")
+# pragma message("!!!!!!!!! CHECK_NNLU set in rspl/rspl.h !!!!!!!!!!!")
#endif
/* Do an arbitrary printf */
@@ -356,7 +359,7 @@ static void rev_test_vram(size_t size) {
static void *rev_malloc(rspl *s, size_t size) {
void *rv;
- if ((size + 1 * 1024 * 1204) > g_test_ram)
+ if ((size + 1 * 1024 * 1024) > g_test_ram)
rev_test_vram(size);
if ((rv = malloc(size)) == NULL) {
rev_reduce_cache(size);
@@ -371,7 +374,7 @@ static void *rev_malloc(rspl *s, size_t size) {
static void *rev_calloc(rspl *s, size_t num, size_t size) {
void *rv;
- if (((num * size) + 1 * 1024 * 1204) > g_test_ram)
+ if (((num * size) + 1 * 1024 * 1024) > g_test_ram)
rev_test_vram(size);
if ((rv = calloc(num, size)) == NULL) {
rev_reduce_cache(num * size);
@@ -386,7 +389,7 @@ static void *rev_calloc(rspl *s, size_t num, size_t size) {
static void *rev_realloc(rspl *s, void *ptr, size_t size) {
void *rv;
- if ((size + 1 * 1024 * 1204) > g_test_ram)
+ if ((size + 1 * 1024 * 1024) > g_test_ram)
rev_test_vram(size);
if ((rv = realloc(ptr, size)) == NULL) {
rev_reduce_cache(size); /* approximation */
@@ -1323,7 +1326,7 @@ set_search_limit(
/* should evaluate in[0..di-1], and return number that is not to exceed */
/* limitv. NULL if not used */
void *lcntx, /* Context passed to limit() */
- double limitv /* Value that limit() is not to exceed */
+ double limitv /* Value that limitf() is not to exceed */
) {
schbase *b = NULL; /* Pointer to search base information structure */
@@ -4025,10 +4028,11 @@ double *err /* Output error (weighted) distance at solution point */
}
#ifdef DEBUG
- if (wsrv2)
+ if (wsrv2) {
DBG(("Got second ink limit triangle in tetrahedron\n"));
- else if (wsrv)
+ } else if (wsrv) {
DBG(("Got first ink limit triangle in tetrahedron\n"));
+ }
#endif
*err = dist;
return wsrv;
@@ -6820,7 +6824,7 @@ static void free_surfhash(rspl *s, int del) {
/* vertex status */
typedef enum {
- vtx_norm = 0, /* Normal vertex in primary bxcell */
+ vtx_norm = 0, /* Normal vertex in primary bxcell - initial value */
vtx_sha = 1, /* Vertex has been shadowed */
vtx_del = 2, /* Vertex has been deleted because it's shadowed */
vtx_oil = 3 /* Vertex is over ink limit */
@@ -6828,8 +6832,9 @@ typedef enum {
struct _vtxrec {
int ix; /* fwd index of vertex */
- int cix; /* Cell index for this vertex */
- double v[MXRO]; /* Output value of vertex */
+ int cix; /* Fwd cell vertex is in index */
+ double vv[MXRO]; /* Output value of vertex */
+ double vl[MXRO]; /* Log compressed output value of vertex */
double dist; /* Distance from center point squared */
int tcount; /* Touch count for converting to fwd cells */
int acount; /* Actual count for converting to fwd cells */
@@ -6866,7 +6871,7 @@ struct _vtxcache {
}; typedef struct _vtxcache vtxcache;
-/* Create the vertex list & hash */
+/* Create the (empty) vertex list & hash */
static void create_vtxrec_list(rspl *s, vtxcache *vc) {
vc->hash_size = primes[3]; /* 3373 */
if ((vc->hash = (vtxrec **) rev_calloc(s, vc->hash_size, sizeof(vtxrec *))) == NULL)
@@ -6949,16 +6954,67 @@ static vtxrec *get_vtxrec(vtxcache *vc, int ix) {
hash = ix % vc->hash_size;
- for (vx = vc->hash[hash]; vx != NULL; vx = vx->hlink) {
+ for (vx = vc->hash[hash]; vx != NULL; vx = vx->hlink) {
if (vx->ix == ix)
return vx;
}
return NULL;
}
+/* Log compress an output value wrt to center point */
+static void logcomp(
+ rspl *s,
+ double *out,
+ double *in,
+ double *cent
+) {
+ int f, fdi = s->fdi;
+ double len;
+
+ if (s->rev.surflin_en) {
+#ifdef NEVER
+ /* (This doesn't seem to improve things) */
+ /* Calculate vector length */
+ for (len = 0.0, f = 0; f < fdi; f++) {
+ double tt= in[f] - cent[f];
+ len += tt * tt;
+ }
+ len = sqrt(len);
+
+ /* change length to log length */
+ if (len > DBL_EPSILON) {
+
+ len = 20.0 * pow(len, 0.25)/len; /* Ratio */
+
+ for (f = 0; f < fdi; f++) {
+ double tt = in[f] - cent[f];
+ out[f] = len * tt + cent[f];
+ }
+ }
+#else
+ if (s->rev.surflin != NULL) {
+ co p;
+
+ for (f = 0; f < fdi; f++)
+ p.p[f] = in[f];
+ s->rev.surflin->interp(s->rev.surflin, &p);
+ for (f = 0; f < fdi; f++)
+ out[f] = p.v[f] - s->rev.linoff[f];
+ } else {
+ for (f = 0; f < fdi; f++)
+ out[f] = in[f];
+ }
+#endif
+ } else {
+ for (f = 0; f < fdi; f++)
+ out[f] = in[f];
+ }
+}
+
/* Create a new vtxrec or return the current one. */
/* Allocates it, adds it to cache. */
/* DOESN"T add it to vtxlist. */
+/* If new, sets status = vtx_norm */
static vtxrec *new_vtxrec(
rspl *s,
vtxcache *vc,
@@ -7005,12 +7061,15 @@ static vtxrec *new_vtxrec(
/* Get the output value */
for (f = 0; f < fdi; f++)
- vx->v[f] = gp[f];
+ vx->vv[f] = gp[f];
+
+ /* Set vl[] */
+ logcomp(s, vx->vl, vx->vv, s->rev.ocent);
/* Compute distance to overall center point squared */
vx->dist = 0.0;
for (f = 0; f < fdi; f++) {
- double tt = gp[f] - s->rev.ocent[f];
+ double tt = vx->vl[f] - s->rev.ocent[f];
vx->dist += tt * tt;
}
@@ -7020,7 +7079,7 @@ static vtxrec *new_vtxrec(
int mi;
double gw = s->rev.gw[f];
double gl = s->rev.gl[f];
- t = (vx->v[f] - gl)/gw;
+ t = (vx->vv[f] - gl)/gw;
mi = (int)floor(t); /* Grid coordinate */
if (mi < 0) /* Limit to valid cube base index range */
mi = 0;
@@ -7442,7 +7501,7 @@ static void extend_bxcell_shadow_group(
/* vertex length from gamut center */
for (vlen= 0.0, f = 0; f < fdi; f++) {
- double tt = vx->v[f] - s->rev.ocent[f];
+ double tt = vx->vl[f] - s->rev.ocent[f];
vlen += tt * tt;
}
vlen = sqrt(vlen);
@@ -7453,7 +7512,7 @@ static void extend_bxcell_shadow_group(
scale = 1.0;
for (f = 0; f < fdi; f++)
- sv[f] = (scale * (vx->v[f] - s->rev.ocent[f])) + s->rev.ocent[f];
+ sv[f] = (scale * (vx->vl[f] - s->rev.ocent[f])) + s->rev.ocent[f];
/* Distance from scaled vertex to group center */
for (w = 0.0, f = 0; f < fdi; f++) {
@@ -8389,13 +8448,13 @@ static void plot_fxcell_surface(rspl *s, int dofclabels, int dobxcells, int dowa
is inside or on the surface of the gamut. A simple and definitive
topological rule hasn't been forthcoming, so a simpler heursitic
of visiblilty from a singe internal "focal" point is currently used.
- calc_ocent() attempts to choose a point with best visibilit of
+ calc_ocent() attempts to choose a point with best visibility of
all the gamut surfaces, since any self-shadowing results in
gamut surface holes.
Improvements would be to create per-axis mappings (and separate
the shadow vertex locations from the real ones) to re-shape
- the gamut into a squere as much as possible.
+ the gamut into a square as much as possible.
Multiple external focal points could be used,
a vertex being shadowed only when it can't be "seen" by any
external focal point. It's hard to figure how to make the latter
@@ -8408,12 +8467,19 @@ static void plot_fxcell_surface(rspl *s, int dofclabels, int dobxcells, int dowa
accelleration structure for shadow testing (BSP tree ??),
and build the gamut surface incrementally from existing
furthest points.
+
+ (Have loop re-orderings been exausted ? i.e. can overlap
+ triangle processing "Do a first pass for each test vertex,
+ testing against just the triangles that are associated with
+ it's triangle" be used for main shadowing testing ?)
*/
/* Struct to slice locus points */
struct _slpoint {
double v[MXRO]; /* Point location */
- double rad; /* Distance from ccent */
+ double wrad; /* Weighted radius */
+ double rad; /* Distance from ccent along slice */
+ double minrad; /* Minimum distance from ccent */
double cvec[MXRO]; /* Vector from this point to ccent */
double len; /* Length of segment, -1 if no good */
@@ -8426,12 +8492,15 @@ struct _ocenctx {
int ares; /* angle resolution */
slpoint *p[MXRO]; /* Slice locus points */
double ccent[MXRO]; /* Construction center point */
+ double ret; /* return value */
+ int oog; /* flag set if center is out of gamut */
int debug;
}; typedef struct _ocenctx ocenctx;
/* Given a set of slice locus points and a proposed center point, */
-/* compute the weighted averag of the orthogonality of the point */
+/* compute the weighted average of the orthogonality of the point */
/* to each locus line segment. (smaller is better) */
+/* (Used for optimizing the focal/center point.) */
static double aorthog(void *_ctx, double *cent) {
ocenctx *ctx = (ocenctx *) _ctx;
rspl *s = ctx->s;
@@ -8441,6 +8510,8 @@ static double aorthog(void *_ctx, double *cent) {
double ang, aang = 0.0;
int naang = 0;
+ ctx->oog = 0;
+
if (ctx->debug) printf("aorthog called with cent %s\n",debPdv(fdi,cent));
for (ff = 0; ff < fdi; ff++) {
@@ -8480,7 +8551,10 @@ static double aorthog(void *_ctx, double *cent) {
/* Normalized difference in distance over length */
/* Compute dot product of cv and segment vector */
+ /* (ang range 0.0 .. 1.0 */
ang = fabs(trad - nrad)/ctx->p[ff][aa].len;
+ if (ang > 1.0)
+ ang = 1.0;
if (ctx->debug) printf(" aa %d: trad %f nrad %f, diff %f, len %f, ang %f\n",aa,trad,nrad,fabs(trad - nrad),ctx->p[ff][aa].len,ang);
@@ -8493,14 +8567,11 @@ static double aorthog(void *_ctx, double *cent) {
if (dot < 0.0) {
if (ctx->debug) printf(" dot is %f\n",dot);
- ang = 50.0; /* Big value */
+ ang = 50.0; /* Big value */
+ ctx->oog = 1;
} else {
ang = pow(ang, 50.0); /* Weight high angles */
}
-
- if (ang > aang)
- aang = ang;
-
aang += ang;
naang++;
}
@@ -8509,6 +8580,8 @@ static double aorthog(void *_ctx, double *cent) {
if (ctx->debug) printf(" returning %f\n",aang);
+ ctx->ret = aang;
+
return aang;
}
@@ -8518,6 +8591,7 @@ static double aorthog(void *_ctx, double *cent) {
/* surface of the gamut. */
static void calc_ocent(rspl *s) {
int i, j, aa, mm;
+ int e, ee, di = s->di;
int f, ff, fdi = s->fdi;
int rgres = s->rev.res; /* number of bwd cells */
double minmax[2][MXRO][MXRO]; /* Range min/max points for each axis */
@@ -8527,14 +8601,15 @@ static void calc_ocent(rspl *s) {
double ss[MXRO];
ocenctx ctx; /* Context */
double atanscale;
+ int ici, nici;
/* Scan the forward array for the min and max points of each axis */
for (f = 0; f < fdi; f++) {
- minmax[0][f][f] = 1e30;
- minmax[1][f][f] = -1e30;
+ minmax[0][f][f] = 1e200;
+ minmax[1][f][f] = -1e200;
}
- /* Scan the Grid for min/max values */
+ /* Scan the fwd Grid for min/max values */
for (gp = s->g.a, ep = s->g.a + s->g.no * s->g.pss; gp < ep; gp += s->g.pss) {
for (ff = 0; ff < fdi; ff++) {
if (minmax[0][ff][ff] > gp[ff]) {
@@ -8550,8 +8625,8 @@ static void calc_ocent(rspl *s) {
if (fdi == 1) {
for (f = 0; f < fdi; f++)
- s->rev.ocent[f] = 0.5 * (minmax[0][0][f] + minmax[1][0][f]);
- DBG(("calc_ocent: got ocent = %s\n",debPdv(fdi,s->rev.ocent)));
+ s->rev.ocent[f] = 0.5 * (minmax[0][f][f] + minmax[1][f][f]);
+ DBG(("calc_ocent: got 1d ocent = %s\n",debPdv(fdi,s->rev.ocent)));
return;
}
@@ -8569,7 +8644,8 @@ static void calc_ocent(rspl *s) {
for (f = 0; f < fdi; f++)
s->rev.ocent[f] = ctx.ccent[f] /= ((fdi-1) * 2.0);
- DBG(("calc_ocent: ccent = %s\n",debPdv(fdi,ctx.ccent)));
+ DBG(("calc_ocent: initial ccent = %s\n",debPdv(fdi,ctx.ccent)));
+//printf("calc_ocent: initial ccent = %s\n",debPdv(fdi,ctx.ccent));
/* If it's all to hard ... */
if (fdi != 3) {
@@ -8581,11 +8657,17 @@ static void calc_ocent(rspl *s) {
midix[f] = (int)((ctx.ccent[f] - s->rev.gl[f])/s->rev.gw[f] + 0.5);
mid[f] = (midix[f]+0.5) * s->rev.gw[f] + s->rev.gl[f];
}
+//printf("calc_ocent: mid point = %s\n",debPdv(fdi,mid));
/* Array for each slice values at angle (+ repeat at end) */
ctx.debug = 0;
ctx.s = s;
- ctx.ares = 20;
+ ctx.ares = (rgres + 1) & ~1; /* Make even so that there is an opposite angle */
+ if (ctx.ares < 6)
+ ctx.ares = 6;
+ else if (ctx.ares > 20)
+ ctx.ares = 20;
+//printf(" ocent ares %d\n",ctx.ares);
atanscale = ctx.ares/(2.0 * DBL_PI);
for (ff = 0; ff < fdi; ff++) {
if ((ctx.p[ff] = (slpoint *)rev_calloc(s, ctx.ares+1,sizeof(slpoint))) == NULL)
@@ -8593,96 +8675,201 @@ static void calc_ocent(rspl *s) {
INCSZ(s, (ctx.ares+1) * sizeof(slpoint));
}
-//printf("~1 locating center point\n");
+ /* Use 5 passes to locate a more reliable initial center point */
+ for (nici = 10, ici = 0; ici < nici; ici++) {
+//printf(" locating center point iter %d\n",ici);
- /* Set initial radius values */
- for (aa = 0; aa < ctx.ares; aa++) {
- for (ff = 0; ff < fdi; ff++)
- ctx.p[ff][aa].rad = -1.0;
- }
-
- /* Take three slices through the rev[] array, plotting */
- /* the maximum circumference for the slice */
+ /* Set initial radius values */
+ for (ff = 0; ff < fdi; ff++) {
+ for (aa = 0; aa < ctx.ares; aa++) {
+ ctx.p[ff][aa].wrad = -1.0;
+ ctx.p[ff][aa].rad = -1.0;
+ ctx.p[ff][aa].minrad = 1e38;
+ }
+ }
+
+ /* Take three slices through the rev[] array, plotting */
+ /* the maximum circumference for the slice */
- /* For the axis we're slicing */
- for (ff = 0; ff < fdi; ff++) {
- DCOUNT(cc, MXRO, 2, 0, 0, rgres); /* Counter through bwd cells */
- double vv[MXRO];
- int aa;
-
+ /* For the axis we're slicing */
+ for (ff = 0; ff < fdi; ff++) {
+ FCOUNT(cc, MXRO, 3); /* Counter through bwd cells */
+ int start[3], endp1[3];
+ double vv[MXRO];
+ int aa;
+
//printf(" slice axis %d\n",ff);
- /* Scan this slice of rev[] */
- DC_INIT(cc);
- while (!DC_DONE(cc)) {
- int ix;
- int slix[2]; /* Indexes in slice direction */
- int *rp;
- int co[MXRO];
-
- /* Compute bx index */
- ix = 0;
- for (j = f = 0; f < fdi; f++) {
- if (f == ff)
- co[f] = midix[f];
- else {
- slix[j] = f;
- co[f] = cc[j++];
+ /* Setup "fat" slice range */
+ for (f = 0; f < fdi; f++) {
+ if (f == ff) {
+ start[f] = midix[f]-1;
+ if (start[f] < 0)
+ start[f] = 0;
+ endp1[f] = midix[f]+2;
+ if (endp1[f] > rgres)
+ endp1[f] = rgres;
+ } else {
+ start[f] = 0;
+ endp1[f] = rgres;
}
- ix += co[f] * s->rev.coi[f];
}
-//printf(" bx %d, %d ix %d\n",co[0],co[1],co[2],ix);
+ FRECONFA(cc, start, endp1);
+
+//printf(" slice range %d - %d, %d - %d, %d - %d\n", start[0], endp1[0]-1, start[1], endp1[1]-1, start[2], endp1[2]-1);
+
+ /* Scan this 3 thick, 2D slice of rev[] */
+ FC_INIT(cc);
+ while (!FC_DONE(cc)) {
+ int ix;
+ int slix[MXRO]; /* Indexes in slice direction + orthogonal */
+ int *rp;
+
+ /* Compute bx index */
+ ix = 0;
+ for (j = f = 0; f < fdi; f++) {
+ ix += cc[f] * s->rev.coi[f];
+ if (f != ff)
+ slix[j++] = f;
+ }
+ slix[j++] = ff;
+//printf(" bx %d, %d ix %d\n",cc[0],cc[1],cc[2],ix);
- if (s->rev.rev[ix] == NULL) {
+ if (s->rev.rev[ix] == NULL) {
//printf(" rev is empty\n");
- goto next_bx;
- }
+ goto next_bx;
+ }
- /* For all the vertex values in bx rev[] */
- for (rp = s->rev.rev[ix]+3; *rp != -1; rp++) {
- float *fcb = s->g.a + *rp * s->g.pss;
- double x, y, rad, ang;
+ /* For all the cubes bx rev[] */
+ for (rp = s->rev.rev[ix]+3; *rp != -1; rp++) {
-//printf(" vtx %d\n",*rp);
- /* Ignore over ink limit values */
- if (s->limiten && fcb[-1] > s->limitv)
- continue;
-
- /* Compute radius and normalize */
- x = fcb[slix[0]] - ctx.ccent[slix[0]];
- y = fcb[slix[1]] - ctx.ccent[slix[1]];
- rad = sqrt(x * x + y * y);
- if (rad < EPS)
- continue;
+ /* For each vertx of this cube */
+ for (ee = 0; ee < (1<<di); ee++) {
+ int vix = *rp + s->g.hi[ee];
+ float *gp = s->g.a + vix * s->g.pss; /* Pointer to float of fwd vertex */
+ double fcb[MXRO];
+ double x, y, z, wrad, rad, ang;
- /* Quantized angle this point is at */
- ang = atanscale * atan2(y, x);
- aa = (int)floor(ang);
- if (aa < 0)
- aa += ctx.ares;
- if (aa >= ctx.ares)
- aa -= ctx.ares;
+ /* Don't add over ink limit vertexes */
+ if (s->limiten && gp[-1] > s->limitv) {
+ continue;
+ }
+
+ for (f = 0; f < fdi; f++)
+ fcb[f] = gp[f];
+
+ /* (Don't) Convert output values to log values */
+ /* logcomp(s, fcb, fcb, ctx.ccent); */
+
+ /* Compute 2D radius and normalize */
+ x = fcb[slix[0]] - ctx.ccent[slix[0]];
+ y = fcb[slix[1]] - ctx.ccent[slix[1]];
+ z = fcb[slix[2]] - ctx.ccent[slix[2]];
+ /* wrad is "elipsoid" weighted radius in slice */
+ wrad = x * x + y * y - 1.5 * z * z;
+ wrad = sqrt(wrad < 0.0 ? 0.0 : wrad);
+ rad = sqrt(x * x + y * y);
+ if (rad < EPS || wrad < EPS)
+ continue;
+
+ /* Quantized angle this point is at */
+ ang = atanscale * atan2(y, x);
+ aa = (int)floor(ang);
+ if (aa < 0)
+ aa += ctx.ares;
+ if (aa >= ctx.ares)
+ aa -= ctx.ares;
//printf(" slice %d vtx %f %f %f rad %f, ang %f aa %d\n", ff, fcb[0], fcb[1], fcb[2], rad, ang, aa);
- if (rad > ctx.p[ff][aa].rad) {
- ctx.p[ff][aa].rad = rad;
+ if (wrad > ctx.p[ff][aa].wrad) {
+ ctx.p[ff][aa].wrad = wrad;
+ ctx.p[ff][aa].rad = rad;
- for (f = 0; f < fdi; f++) {
- ctx.p[ff][aa].v[f] = fcb[f];
+ /* Copy far point */
+ for (f = 0; f < fdi; f++)
+ ctx.p[ff][aa].v[f] = fcb[f];
+
+ /* (don't) Flatten the points to lie on the notional center */
+ /* ctx.p[ff][aa].v[ff] = ctx.ccent[ff]; */
+ }
+
+ /* Track min in case ccent is not within slice */
+ if (rad < ctx.p[ff][aa].minrad) {
+ ctx.p[ff][aa].minrad = rad;
+ }
}
- /* Flatten the points to lie on the notional center */
- ctx.p[ff][aa].v[ff] = ctx.ccent[ff];
}
+ next_bx:;
+ FC_INC(cc);
}
- next_bx:;
- DC_INC(cc);
+
+ /* Repeat first in extra at end */
+ ctx.p[ff][ctx.ares] = ctx.p[ff][0]; /* Structure copy */
}
- /* Repeat first in extra at end */
- ctx.p[ff][ctx.ares] = ctx.p[ff][0]; /* Structure copy */
+ /* Check if center point is within slice by looking for empty entries. */
+ {
+ double ccvec[MXRO]; /* Center correction vector */
+ double ccount = 0.0;
+
+ for (f = 0; f < fdi; f++)
+ ccvec[f] = 0.0;
+
+ for (ff = 0; ff < fdi; ff++) {
+ for (aa = 0; aa < ctx.ares; aa++) {
+
+//printf(" slice %d aa %d, vtx %s rad %f, irad %f\n", ff, aa, debPdv(fdi,ctx.p[ff][aa].v), ctx.p[ff][aa].rad, ctx.p[ff][aa].minrad);
+
+ /* Either the grid is very sparse, or our center */
+ /* is outside */
+ if (ctx.p[ff][aa].rad < 0.0) {
+ int oaa = aa + (ctx.ares/2);
+ if (oaa >= ctx.ares)
+ oaa -= ctx.ares;
+
+//printf(" oaa %d, vtx %s rad %f, irad %f\n", oaa, debPdv(fdi,ctx.p[ff][oaa].v), ctx.p[ff][oaa].rad, ctx.p[ff][oaa].minrad);
+
+ /* If oposite side has an entry */
+ if (ctx.p[ff][oaa].rad > 0.0) {
+ double cor[MXRO];
+ double prop = (3.0 * ctx.p[ff][oaa].minrad
+ + ctx.p[ff][oaa].rad)/4.0;
+
+ prop /= ctx.p[ff][oaa].rad; /* Proportion of distance to v[] */
+
+ for (f = 0; f < fdi; f++)
+ cor[f] = prop * (ctx.p[ff][oaa].v[f] - ctx.ccent[f]);
+//printf(" prop %f, corr %s\n",prop,debPdv(fdi,cor));
+ for (f = 0; f < fdi; f++)
+ ccvec[f] += cor[f];
+ ccount++;
+ }
+ }
+ }
+ }
+
+//printf("ccount %f\n",ccount);
+ if (ccount > 0.0) { /* Make adjustment */
+ if (ici < (nici-1)) {
+ for (f = 0; f < fdi; f++)
+ ccvec[f] /= ccount;
+//printf("Corecting center by %s\n",debPdv(fdi,ccvec));
+ for (f = 0; f < fdi; f++)
+ ctx.ccent[f] += ccvec[f];
+//printf("cceny now %s\n",debPdv(fdi,ctx.ccent));
+ } else { /* Last round and correction needed */
+ if (0.0 && s->verbose)
+ fprintf(stdout, "%cFailed to locate aprox. gamut center\n",cr_char);
+ }
+ } else {
+ break; /* We're done */
+ }
+ }
}
+//printf("calc_ocent: refined ccent = %s\n",debPdv(fdi,ctx.ccent));
+
/* Pre-compute point to point info to speed optimization */
for (ff = 0; ff < fdi; ff++) {
for (aa = 0; aa < ctx.ares; aa++) {
@@ -8706,7 +8893,9 @@ static void calc_ocent(rspl *s) {
/* slice segment. This should maximize visibility of the inner of the */
/* gamut surface, for shadow testing. */
for (f = 0; f < fdi; f++)
- ss[f] = 5.0;
+ ss[f] = fabs(0.1 * (minmax[1][f][f] - minmax[0][f][f]));
+
+//ctx.debug = 1;
/* return 0 on sucess, 1 on failure due to excessive itterations */
if (powell(NULL, fdi, s->rev.ocent, ss, 1e-3, 500, aorthog, (void *)&ctx, NULL, NULL)) {
@@ -8715,8 +8904,16 @@ static void calc_ocent(rspl *s) {
s->rev.ocent[f] = ctx.ccent[f];
}
-// ctx.debug = 1;
-// printf("Final angle = %f\n", aorthog(&ctx, ctx.ccent));
+//ctx.debug = 1;
+
+ /* Check result */
+ aorthog(&ctx, ctx.ccent);
+
+ /* Hmm. This isn't very reliable in detecting failure. */
+ if (ctx.oog)
+ printf("calc_ocent failed to return in-gamut focal point!\n");
+
+//printf("Final angle = %f\n", ctx.ret);
#ifdef REVVRML /* Plotting routine declarations */
/* Diagnostic - dump the gamut slice locii */
@@ -8729,7 +8926,7 @@ static void calc_ocent(rspl *s) {
double blue[3] = { 0.1, 0.1, 0.8 };
double *rgb[3] = { red, green, blue };
- wrl = new_vrml("section", 0, vrml_lab);
+ wrl = new_vrml("section", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
/* Show vertex labels */
@@ -8772,7 +8969,384 @@ static void calc_ocent(rspl *s) {
DECSZ(s, (ctx.ares+1) * sizeof(slpoint));
}
- DBG(("calc_ocent: got ocent = %s\n",debPdv(fdi,s->rev.ocent)));
+ DBG(("calc_ocent: final ocent = %s\n",debPdv(fdi,s->rev.ocent)));
+//printf("calc_ocent: final ocent = %s\n",debPdv(fdi,s->rev.ocent));
+}
+
+/* Create gamut surface linearization (surflin) transform. */
+/* This is used by logcomp() to try and straighten out the */
+/* device response so that the ocent is "visible" from */
+/* any point on the surface. */
+/* (We assume we are called at the correct point when bx->status == bx_uninit) */
+static int calc_surflin(
+rspl *s,
+vtxcache *vc, /* Vertexes */
+assdire *edgdir /* Edge lookup for vertex */
+) {
+ int i, j, g;
+ int e, ee, di = s->di;
+ int f, ff, fdi = s->fdi;
+ vtxrec *vx, *nvx;
+ int nitter, itter;
+
+ double vxv[POW2MXRI][MXRO]; /* Overal fwd interp vertex values */
+ int nvtx;
+ cow *mpoints;
+ int gres[MXDO];
+ double min[MXDO], max[MXDO];
+ double vmin[MXDO], vmax[MXDO];
+
+ for (f = 0; f < fdi; f++)
+ gres[f] = s->g.bres; // ??
+
+//printf("calc_surflin: rspl res %d\n",gres[0]);
+ DBG(("calc_surflin: rspl res %d\n",gres[0]));
+
+//printf("~1 gres = %d %d %d\n", s->g.res[0], s->g.res[1], s->g.res[2]);
+//printf("~1 ci = %d %d %d\n", s->g.ci[0], s->g.ci[1], s->g.ci[2]);
+
+ /* Lookup interpolation cube corners */
+ for (ee = 0; ee < (1<<di); ee++) {
+ int ix;
+ float *fcb;
+ for (ix = e = 0; e < di; e++) {
+ if (ee & (1<<e))
+ ix += s->g.ci[e] * (s->g.res[e] - 1);
+ }
+ fcb = s->g.a + ix * s->g.pss;
+ for (f = 0; f < fdi; f++)
+ vxv[ee][f] = fcb[f];
+//printf("~1 cube corners %d = %f %f %f\n",ee, vxv[ee][0], vxv[ee][1], vxv[ee][2]);
+ }
+
+//printf("calc_surflin: counting number of vertexes:\n");
+ /* Count the number of vertexes we may need */
+ nvtx = 0;
+ for (i = 0; i < vc->hash_size; i++) {
+ for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+ nvtx++;
+ }
+ }
+ DBG(("calc_surflin: %d mapping points\n",nvtx));
+//printf("calc_surflin: %d mapping points\n",nvtx);
+
+//printf("calc_surflin: computing goal values:\n");
+ /* Currently the vertex vl = vv = output value of vertex. */
+ /* Temporarily replace vv with the idealize (linear interp) output "goal" values. */
+ for (i = 0; i < vc->hash_size; i++) {
+ for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+ int tix; /* Temp fwd cell index */
+ double we[MXRI]; /* Vertex input position == 1.0 - Weight */
+ double gw[POW2MXRI]; /* weight for each grid cube corner */
+ double w;
+
+ /* Compute this vertexes relative input position */
+ for (tix = vx->ix, e = 0; e < di; e++) {
+ int dix;
+ dix = tix % s->g.res[e];
+ tix /= s->g.res[e];
+ we[e] = (double)dix/(s->g.res[e]-1.0);
+ }
+
+ /* Compute corner weights needed for interpolation */
+ gw[0] = 1.0;
+ for (e = 0, g = 1; e < di; g *= 2, e++) {
+ for (j = 0; j < g; j++) {
+ gw[g+j] = gw[j] * we[e];
+ gw[j] *= (1.0 - we[e]);
+ }
+ }
+
+ /* Linear interpolated output values */
+ w = gw[0];
+ for (f = 0; f < fdi; f++) /* Base of cube */
+ vx->vv[f] = w * vxv[0][f];
+
+ for (g = 1; g < (1<<di); g++) { /* For all other corners of cube */
+ w = gw[g];
+ for (f = 0; f < fdi; f++)
+ vx->vv[f] += w * vxv[g][f];
+ }
+ }
+ }
+
+ /* Now itteratively adjust the vl values to better match the scaled */
+ /* relative positions of the goal values. */
+
+ /* Go through them again to get every line they are part of */
+ /* (We're assuming we need exponentially more itters with finer point */
+ /* spacing ?) */
+ nitter = (int)(0.06 * s->g.bres * s->g.bres + 0.5);
+ if (nitter < 1)
+ nitter = 1;
+ for (itter = 0; itter < nitter; itter++) {
+
+ DBG(("calc_surflin: maping itter %d\n",itter));
+//printf("calc_surflin: maping itter %d\n",itter);
+
+ for (i = 0; i < vc->hash_size; i++) {
+ for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+ assdire *edg; /* Edge table */
+ float *fp;
+ int fl;
+ double agrad, aorad;
+ int nn;
+ double scale;
+
+ fp = s->g.a + vx->ix * s->g.pss; /* This vertex in fwd grid */
+ fl = FLV(fp); /* Edge flags for this vertex */
+ edg = edgdir + fl;
+
+ /* For vertexes at the end of all possible edges common with this vertex, */
+ /* compute average radius */
+ agrad = aorad = 0.0;
+ for (j = 0; j < edgdir[fl].no; j++) {
+ int eix;
+
+ /* Index number of vertex other than the one we got it from */
+ if (edg->ti[j].goffs[0] != 0)
+ eix = vx->ix + edg->ti[j].goffs[0];
+ else
+ eix = vx->ix + edg->ti[j].goffs[1];
+
+ if ((nvx = get_vtxrec(vc, eix)) != NULL) {
+ double glen, olen;
+ glen = olen = 0.0;
+ for (f = 0; f < fdi; f++) {
+ double tt;
+
+ tt = vx->vv[f] - nvx->vv[f];
+ glen += tt * tt;
+
+ tt = vx->vl[f] - nvx->vl[f];
+ olen += tt * tt;
+ }
+ glen = sqrt(glen);
+ olen = sqrt(olen);
+ agrad += glen;
+ aorad += olen;
+ nn++;
+ }
+ }
+
+ if (nn == 0) { /* Hmm. No neighors ? */
+ vx->status = vtx_del; /* Mark it as isolated */
+ continue;
+ }
+
+ scale = aorad/agrad; /* Local scale factor from goal to output */
+
+ /* Reset the current vertex output value based on the relative */
+ /* position of it in goal space */
+ for (f = 0; f < fdi; f++)
+ vx->vl[f] = 0.0;
+
+ nn = 0;
+ for (j = 0; j < edgdir[fl].no; j++) {
+ int eix;
+
+ /* Index number of vertex other than the one we got it from */
+ if (edg->ti[j].goffs[0] != 0)
+ eix = vx->ix + edg->ti[j].goffs[0];
+ else
+ eix = vx->ix + edg->ti[j].goffs[1];
+
+ if ((nvx = get_vtxrec(vc, eix)) != NULL) {
+ for (f = 0; f < fdi; f++)
+ vx->vl[f] += nvx->vl[f] + scale * (vx->vv[f] - nvx->vv[f]);
+ nn++;
+ }
+ }
+ for (f = 0; f < fdi; f++)
+ vx->vl[f] /= (double)nn;
+ }
+ }
+ }
+
+ DBG(("calc_surflin: creating rspl\n"));
+//printf("calc_surflin: creating rspl\n");
+
+ /* Now construct rspl setup mapping points from vertex normal output values */
+ /* to adjusted vl values */
+
+ if ((s->rev.surflin = new_rspl(RSPL_NOFLAGS, fdi, fdi)) == NULL)
+ error("calc_surflin: new_rspl failed");
+
+ nvtx++; /* One for center point */
+
+ /* Allocate rspl setup points */
+ if ((mpoints = malloc(sizeof(cow) * 2 * nvtx)) == NULL)
+// if ((mpoints = malloc(sizeof(cow) * nvtx)) == NULL)
+ error("calc_surflin: malloc of %d rspl setup points failed",nvtx);
+
+ nvtx = 0;
+
+#ifndef NEVER
+ /* Center point */
+ for (f = 0; f < fdi; f++) {
+ mpoints[nvtx].p[f] = s->rev.ocent[f];
+ mpoints[nvtx].v[f] = s->rev.ocent[f];
+ }
+ mpoints[nvtx].w = 10.0;
+ nvtx++;
+#endif
+
+ /* Set the surface mapping points and restore vertexes values */
+ for (i = 0; i < vc->hash_size; i++) {
+ for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+ float *fcb;
+
+ fcb = s->g.a + vx->ix * s->g.pss; /* This vertex in fwd grid */
+
+ /* Actual output value as source of mapping */
+ for (f = 0; f < fdi; f++)
+ mpoints[nvtx].p[f] = vx->vv[f] = fcb[f];
+
+ if (vx->status != vtx_norm) { /* Skip isolated values */
+ vx->status = vtx_norm; /* Restore vtx contents */
+ for (f = 0; f < fdi; f++)
+ vx->vl[f] = vx->vv[f];
+ continue;
+ }
+
+ for (f = 0; f < fdi; f++)
+ mpoints[nvtx].v[f] = vx->vl[f];
+ mpoints[nvtx].w = 1.0;
+ nvtx++;
+
+ vx->status = vtx_norm; /* Restore vtx contents */
+ for (f = 0; f < fdi; f++)
+ vx->vl[f] = vx->vv[f];
+
+//printf("Map[%d] %f %f %f -> %f %f %f\n", nvtx-1, mpoints[nvtx-1].p[0], mpoints[nvtx-1].p[1], mpoints[nvtx-1].p[2], mpoints[nvtx-1].v[0], mpoints[nvtx-1].v[1], mpoints[nvtx-1].v[2]);
+
+#ifndef NEVER
+ /* Add intermediate "fixed" point */
+ for (f = 0; f < fdi; f++) {
+ mpoints[nvtx].p[f] = mpoints[nvtx].v[f]
+ = 0.5 * mpoints[nvtx-1].p[f] + 0.5 * s->rev.ocent[f];
+ }
+ mpoints[nvtx].w = 0.5;
+ nvtx++;
+#endif
+ }
+ }
+
+ for (f = 0; f < fdi; f++) {
+ min[f] = 1e38;
+ max[f] = -1e38;
+ vmin[f] = 1e38;
+ vmax[f] = -1e38;
+ }
+ for (i = 0; i < nvtx; i++) {
+
+#ifdef NEVER
+ /* Blend with original values */
+ for (f = 0; f < fdi; f++)
+ mpoints[i].v[f] = 0.5 * mpoints[i].p[f] + 0.5 * mpoints[i].v[f];
+#endif
+
+ for (f = 0; f < fdi; f++) {
+ if (mpoints[i].p[f] < min[f])
+ min[f] = mpoints[i].p[f];
+ if (mpoints[i].p[f] > max[f])
+ max[f] = mpoints[i].p[f];
+
+ if (mpoints[i].v[f] < vmin[f])
+ vmin[f] = mpoints[i].v[f];
+ if (mpoints[i].v[f] > vmax[f])
+ vmax[f] = mpoints[i].v[f];
+ }
+ }
+
+#ifdef REVVRML /* Plot mapping vectors red->green */
+ {
+ vrml *wrl;
+ double red[3] = { 1.0, 0.0, 0.0 };
+ double green[3] = { 0.0, 1.0, 0.0 };
+
+ wrl = new_vrml("suflinvecss", 0, vrml_lab);
+ wrl->start_line_set(wrl, 0);
+
+ for (i = 0; i < nvtx; i++) {
+ wrl->add_col_vertex(wrl, 0, mpoints[i].p, red);
+ wrl->add_col_vertex(wrl, 0, mpoints[i].v, green);
+
+ }
+
+ wrl->make_lines(wrl, 0, 2);
+ wrl->del(wrl);
+ }
+#endif
+
+ DBG(("calc_surflin: mapping points set, about to creat rspl:\n"));
+//printf("calc_surflin: mapping points set, about to creat rspl:\n");
+
+ /* Fit the rspl */
+ s->rev.surflin->fit_rspl_w(s->rev.surflin, RSPL_NOFLAGS, mpoints, nvtx,
+ min, max, gres, vmin, vmax, 4.0, NULL, NULL);
+
+ DBG(("calc_surflin: mapping created\n"));
+//printf("calc_surflin: mapping created\n");
+
+#ifdef NEVER
+ {
+ double de;
+ co p;
+ extern double icmNorm33(double *, double *);
+
+ /* Check fit */
+ de = 0.0;
+
+ for (i = 0; i < nvtx; i++) {
+ for (f = 0; f < fdi; f++)
+ p.p[f] = mpoints[i].p[f];
+
+ s->rev.surflin->interp(s->rev.surflin, &p);
+
+ de += icmNorm33(mpoints[i].v, p.v);
+ }
+ de = de/(double)nvtx;
+
+ printf("Avg fit error = %f\n",de);
+ }
+#endif
+
+ free(mpoints);
+
+ /* Lookup ocent mapping offset */
+ {
+ co p;
+
+ for (f = 0; f < fdi; f++)
+ p.p[f] = s->rev.ocent[f];
+ s->rev.surflin->interp(s->rev.surflin, &p);
+//printf("opoint mapping %f %f %f -> %f %f %f\n", p.p[0], p.p[1], p.p[2], p.v[0], p.v[1], p.v[2]);
+ for (f = 0; f < fdi; f++)
+ s->rev.linoff[f] = p.v[f] - s->rev.ocent[f];
+ }
+
+
+#ifndef NEVER
+ /* Put the transform into use */
+ s->rev.surflin_en = 1;
+
+ /* Transform all the vertexes */
+ for (i = 0; i < vc->hash_size; i++) {
+ for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+ logcomp(s, vx->vl, vx->vv, s->rev.ocent);
+
+ /* Compute distance to overall center point squared */
+ vx->dist = 0.0;
+ for (f = 0; f < fdi; f++) {
+ double tt = vx->vl[f] - s->rev.ocent[f];
+ vx->dist += tt * tt;
+ }
+ }
+ }
+#endif
+
+ return 0;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
@@ -8800,6 +9374,7 @@ rspl *s
/* Note that bit 1 can be set for cells that are not */
/* to be explored because they are in the gamut interior, */
/* and because they have already been added to the seedlist. */
+ int pass; /* Construction pass */
float *gp; /* Pointer to fwd grid points */
DCOUNT(gg, MXRO, fdi, 0, 0, rgres); /* Track the prime seed coordinate */
@@ -8835,9 +9410,24 @@ rspl *s
int nnmxrevcelldepth = 0; /* Maximum nnrev[] list lengths */
int nnrevshare = 0; /* Sum of nnrev[] list reference counts */
#endif
+ datao rgmin, rgmax;
DBG(("init_revaccell called, di = %d, fdi = %d, mgres = %d\n",di,fdi,(int)s->g.mres));
+ /* To help VRML diagnostics, make a guess as to whether the output */
+ /* space is XYZ like, or L*a*b* like */
+
+ s->get_out_range(s, rgmin, rgmax); /* overall output min/max */
+
+ if (fdi >= 3
+ && rgmin[0] >= -1.0 && rgmax[0] < 3.0
+ && rgmin[1] >= -1.0 && rgmax[1] < 3.0
+ && rgmin[2] >= -1.0 && rgmax[2] < 3.0) {
+ s->rev.probxyz = 1;
+ if (s->verbose)
+ fprintf(stdout, "%cLooks like an XYZ space\n",cr_char);
+ }
+
if (fdi > 1 && s->verbose)
fprintf(stdout, "%cInitializing nnrev arrays...\n",cr_char);
@@ -8885,7 +9475,7 @@ rspl *s
/*
* The rev[] and nnrev[] grids contain pointers to lists of grid cube base indexes.
* If the pointer is NULL, then there are no base indexes in that list.
- * A non NULL list uses element [0] to indicate the alocation size of the list,
+ * A non NULL list uses element [0] to indicate the allocation size of the list,
* [1] contains the index of the next free location, [2] contains the reference
* count (lists may be shared), the list starts at [3]. The last entry is marked with -1.
*/
@@ -8897,6 +9487,8 @@ rspl *s
double iv[MXDI]; /* Input value corresponding to grid */
DBG(("Looking up fwd vertex ink limit values\n"));
+//printf("Looking up fwd vertex ink limit values\n");
+//printf("s->limitv = %f\n",s->limitv);
/* Calling the limit function for each fwd vertex could be bad */
/* if the limit function is slow. Maybe an octree type algorithm */
/* could be used if this is a problem ? */
@@ -8906,7 +9498,9 @@ rspl *s
for (e = 0; e < di; e++)
iv[e] = s->g.l[e] + gc[e] * s->g.w[e]; /* Input sample values */
gp[-1] = (float)(INKSCALE * s->limitf(s->lcntx, iv));
+//printf("~1 set ix %d limitv to %f\n",i,gp[-1]);
}
+//else printf("~1 ix %d limitv is %f\n",i,gp[-1]);
EC_INC(gc);
}
s->g.limitv_cached = 1;
@@ -8958,6 +9552,7 @@ rspl *s
uil = oil = 0;
for (f = 0; f < fdi; f++) /* Init output min/max */
min[f] = max[f] = gp[f];
+
if (!s->limiten || gp[-1] <= s->limitv)
uil = 1;
else
@@ -9065,7 +9660,7 @@ rspl *s
} /* Next reverse grid point in intersecting cube */
} /* Next base grid point */
- DBG(("We skipped %d cells that were over the limit\n",nskcells));
+ DBG(("We skipped %d/%d cells that were over the limit\n",nskcells,gno));
#ifdef CHECK_NNLU
if (fdi > 1) {
@@ -9180,7 +9775,6 @@ rspl *s
cc[f] += (ff & 1) ? 1 : -1;
vflagp += (ff & 1) ? s->rev.coi[f] : -s->rev.coi[f];
-
/* Out of bounds or empty */
if (cc[f] < 0 || cc[f] >= rgres || ((*vflagp & 0xf) == 0)) {
vflag[i] = (vflag[i] & ~0xf) | 2; /* Convert this one to empty surface cell */
@@ -9256,7 +9850,7 @@ rspl *s
}
}
- } else {
+ } else { /* di >= 3 gamut surface finding */
#ifdef REVVRML
/* Plot the initial surface bxcells & their fwd cells. */
@@ -9289,7 +9883,7 @@ rspl *s
/* - - - - - - - - - - - - - - */
/* fill, thin and add, until */
/* there is no more work to do. */
- for (;;) {
+ for (pass = 0;; pass++) {
int phase;
int morevtxadded = 0;
#if defined(REVTABLESTATS) || defined(DEBUG)
@@ -9301,7 +9895,7 @@ rspl *s
/* For each surface bxcell, convert the corresponding */
/* rev[] fwd cubes into vertices. */
- /* (Musk keep bxcells even if none of their verticies */
+ /* (Must keep bxcells even if none of their verticies */
/* are physically in them, so that those verticies get thinned. */
/* could only remove them if vertex was not in any surface cell ?) */
for (pbx = &s->rev.surflist, bx = *pbx; bx != NULL; bx = nbx) {
@@ -9368,10 +9962,10 @@ rspl *s
/* Expand a bxcell's shadow testing group values based on it's vertex list */
/* so that shadow testing works correctly for vertexes that don't */
/* actually lie within the bxcell. (Note that in fact the triangle */
- /* testing creates triangles that are mode of vertexes that may not */
- /* be in this bx's list, so the shadow size doesn't accturatly reprsent */
+ /* testing creates triangles that are made of vertexes that may not */
+ /* be in this bx's list, so the shadow size doesn't accuratly reprsent */
/* the possible shadow area. It's not clear what consequences this has, */
- /* if any. If we extanded the group to cover this, we would need to have ) */
+ /* if any. If we extanded the group to cover this, we would need to have */
/* two groups, a shadower group including those vertexes, and a shadowee */
/* goup for just those vertexes that are part of the bx. */
extend_bxcell_shadow_group(s, &vc, bx);
@@ -9380,6 +9974,19 @@ rspl *s
morevtxadded = 1;
}
+ /* Compute transform rspl that helps "unfold" any regions of the surface */
+ /* that overlap from the perspective of ocent, to try and avoid gaps in */
+ /* the final gamut surface. Existing vtxrec are converted to have vl */
+ /* in the unfolded space. */
+ if (pass == 0 /* && function flag set */) {
+#ifndef EN_UNTWIST /* Control using an environment variable */
+ if (getenv("ARGYLL_UNTWIST_GAMUT_SURFACE") != NULL)
+#endif
+ {
+ calc_surflin(s, &vc, edgdir);
+ }
+ }
+
DBG(("thinning surface vertex lists and converting to cells\n"));
/* (Sorting bxcells doesn't seem to make any performace difference.) */
@@ -9440,8 +10047,8 @@ rspl *s
/* If any of bx is further from nbx and their bounding */
/* cylinders overlap in perspective from rev.ocenter, */
/* assume nbx is a shadow */
- if (shadow_group_group(s, s->rev.ocent, bx->g.bcent, bx->cc, bx->dw,
- nbx->g.bcent, nbx->cc, nbx->dw)) {
+ if (shadow_group_group(s, s->rev.ocent, bx->g.bcent, bx->cc,
+ bx->dw, nbx->g.bcent, nbx->cc, nbx->dw)) {
nbx->wlist = bx->wlist;
bx->wlist = nbx;
//printf("~1 adding shadow nnrev[%d] from xlist\n",nbx->ix);
@@ -9511,8 +10118,8 @@ rspl *s
/* marked deleted ??) */
if (
// vx->status == vtx_norm &&
- shadow_group_vertex(s, s->rev.ocent, bx->g.bcent, bx->cc, bx->dw,
- vx->v)) {
+ shadow_group_vertex(s, s->rev.ocent, bx->g.bcent, bx->cc,
+ bx->dw, vx->vl)) {
add_vtxrec_list(&vc, vx, 0); /* Add if not deleted */
//printf(" Added ix %d from bx %d\n",vx->ix,nbx->ix);
}
@@ -9529,7 +10136,8 @@ rspl *s
error("Failed to find vertex %s in cache",*rp);
if (vx->status == vtx_norm &&
- shadow_group_vertex(s, s->rev.ocent, bx->g.bcent, bx->cc, bx->dw, vx->v)) {
+ shadow_group_vertex(s, s->rev.ocent, bx->g.bcent, bx->cc,
+ bx->dw, vx->vl)) {
add_vtxrec_list(&vc, vx, 1); /* Add if not hidden/deleted */
}
}
@@ -9543,7 +10151,7 @@ rspl *s
/* For vertexes of this bxcell and shadowers, */
/* in order from largst to smallest distance from center. */
for (vx = vc.vtxlist; vx != NULL; vx = vx->tlist) {
- float *vp; /* Vertex being tested */
+ float *fcb; /* Vertex being tested */
int fl;
assdire *tri; /* Triangle table */
@@ -9556,8 +10164,8 @@ rspl *s
//printf("~1 doing vertex %d at %s dist %f\n",vx->ix, debPdv(fdi,vx->v), sqrt(vx->dist));
- vp = s->g.a + vx->ix * s->g.pss; /* This vertex in fwd grid */
- fl = FLV(vp); /* Edge flags for this vertex */
+ fcb = s->g.a + vx->ix * s->g.pss; /* This vertex in fwd grid */
+ fl = FLV(fcb); /* Edge flags for this vertex */
tri = tridir + fl;
//printf("~1 fl %d = 0o%o, no triangles %d\n",fl, fl, tri->no);
@@ -9597,7 +10205,7 @@ rspl *s
for (e = 0; e <= sdi; e++) {
for (f = 0; f < fdi; f++)
- v[e][f] = trivx[e]->v[f];
+ v[e][f] = trivx[e]->vl[f];
}
/* Compute shadow group params of triangle for quick vertex test */
@@ -9633,7 +10241,7 @@ rspl *s
/* Do quick check against triangle */
if (!shadow_group_vertex(s,
- s->rev.ocent, gc, cc, dw, nvx->v)) {
+ s->rev.ocent, gc, cc, dw, nvx->vl)) {
//printf("~1 shadow group check shows no intersection\n");
continue;
}
@@ -9644,11 +10252,14 @@ rspl *s
shdwd = wsrv = 0;
/* Compute line delta */
- vp = s->g.a + nvx->ix * s->g.pss;
- for (f = 0; f < fdi; f++) {
- pv[f] = vp[f];
+ fcb = s->g.a + nvx->ix * s->g.pss;
+ for (f = 0; f < fdi; f++)
+ pv[f] = fcb[f];
+
+ logcomp(s, pv, pv, s->rev.ocent);
+
+ for (f = 0; f < fdi; f++)
de[f] = pv[f] - s->rev.ocent[f];
- }
/* Setup line cla and clb */
init_line_eq_imp(s, NULL, &cla, clb, s->rev.ocent, de, 0);
@@ -9773,13 +10384,13 @@ rspl *s
/* Check all of its neighbor vertexes, to see if */
/* it's safe to actually delete them. */
if (vx->status >= vtx_sha) { /* vertex to delete ? */
- float *vp;
+ float *fcb;
int fl;
assdire *edg; /* Edge table */
//printf("Checking vx %d neighbors\n",vx->ix);
- vp = s->g.a + vx->ix * s->g.pss; /* This vertex in fwd grid */
- fl = FLV(vp); /* Edge flags for this vertex */
+ fcb = s->g.a + vx->ix * s->g.pss; /* This vertex in fwd grid */
+ fl = FLV(fcb); /* Edge flags for this vertex */
edg = edgdir + fl;
/* For all possible edges that use this vertex */
@@ -9972,7 +10583,7 @@ rspl *s
/* Struct to hold test vertex locations */
struct _tvxrec {
- double v[MXRO];
+ double v[MXRO]; /* Log output vertex value */
double dist; /* Distance from center point squared */
int ix[MXRO+1]; /* Indexes of the triangle verticies */
int shad; /* Test result */
@@ -10082,10 +10693,10 @@ rspl *s
for (j = 0; j <= sdi; j++) {
if (trivx[j]->status == vtx_norm) {
for (f = 0; f < fdi; f++)
- tvx->v[f] += 0.95/nntvsh * trivx[j]->v[f];
+ tvx->v[f] += 0.95/nntvsh * trivx[j]->vl[f];
} else {
for (f = 0; f < fdi; f++)
- tvx->v[f] += 0.05/ntvsh * trivx[j]->v[f];
+ tvx->v[f] += 0.05/ntvsh * trivx[j]->vl[f];
trivx[j]->cross = 1; /* For diagnostics */
}
}
@@ -10202,7 +10813,7 @@ rspl *s
for (j = 0; j <= sdi; j++) {
for (f = 0; f < fdi; f++)
- v[j][f] = trivx[j]->v[f];
+ v[j][f] = trivx[j]->vl[f];
}
/* Compute shadow group params of triangle for quick vertex test */
@@ -10380,7 +10991,7 @@ rspl *s
for (j = 0; j <= sdi; j++) {
for (f = 0; f < fdi; f++)
- v[j][f] = trivx[j]->v[f];
+ v[j][f] = trivx[j]->vl[f];
}
/* Compute shadow group params of triangle for quick vertex test */
@@ -10570,7 +11181,7 @@ rspl *s
vtxrec *vx;
if ((vx = get_vtxrec(&vc, *rp)) == NULL)
- continue; /* Hmm. */
+ continue; /* Hmm. Delete it.*/
/* Keep all the un-shadowed or preserved vertexes */
if (vx->status == vtx_norm
@@ -10618,9 +11229,9 @@ rspl *s
/* Add extra over ink limit vertexes. */
for (bx = s->rev.surflist; bx != NULL; bx = bx->slist) {
int sdi = 1; /* sub-simplexes are edges */
- int *crp, *rp, *nrp;
- int ttouch;
+ int *rp;
vtxrec *vx, *nvx;
+ int *exlist = NULL;
/* Add over ink limit vertexes, so that fwd cells will straddle */
/* the ink limit boundary. */
@@ -10630,8 +11241,8 @@ rspl *s
/* over ink limit verticies.) */
if (s->limiten && vflag[bx->ix] & 0x10) {
int *rp;
-//printf("~1 ink limitin is enabled bx %d\n", bx->ix);
+//printf("~1 ink limitin is enabled bx %d\n", bx->ix);
for (rp = bx->sl+3; *rp != -1; rp++) {
float *vp, *evp;
int fl;
@@ -10671,7 +11282,7 @@ rspl *s
//printf(" Checking edge %d (%f) -> %d (%f)\n", vx->ix, vp[-1], eix, evp[-1]);
- /* If over limit, add it to the list */
+ /* If over limit, add it to the expansion list */
if (evp[-1] > s->limitv) {
//printf("~1 added over ink limit vertex %d\n",eix);
@@ -10679,10 +11290,18 @@ rspl *s
continue; /* Added by another bx */
nvx = new_vtxrec(s, &vc, eix);
nvx->status = vtx_oil;
- add2indexlist(s, &bx->sl, eix, 0);
+ add2indexlist(s, &exlist, eix, 0);
}
}
}
+
+ /* If we found over ink limit verticies, add them to our list */
+ if (exlist != NULL) {
+ for (rp = exlist+3; *rp != -1; rp++) {
+ add2indexlist(s, &bx->sl, *rp, 0);
+ }
+ free_indexlist(s, &exlist);
+ }
}
}
@@ -10808,6 +11427,11 @@ rspl *s
bx->status = bx_conv;
}
+ if (s->rev.surflin != NULL) { /* Don't need surflin anymore */
+ s->rev.surflin->del(s->rev.surflin);
+ s->rev.surflin = NULL;
+ s->rev.surflin_en = 0;
+ }
if (cla != NULL)
free_dmatrix(cla, 0, fdi-1, 0, fdi);
free_trirec(s, &stc);
@@ -11911,7 +12535,7 @@ int dofwlabels /* Plot fwd cell base indexs */
double grey[3] = { 0.5, 0.5, 0.5 };
double white[3] = { 1.0, 1.0, 1.0 };
- wrl = new_vrml("raw_bxfwcells", 0, vrml_lab);
+ wrl = new_vrml("raw_bxfwcells", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
if (dofwlabels) {
@@ -12080,7 +12704,7 @@ double xv[MXRO] /* Intersection point */
double blue[3] = { 0.1, 0.1, 0.8 };
double yellow[3] = { 0.8, 0.8, 0.1 };
- wrl = new_vrml("tri_check", 0, vrml_lab);
+ wrl = new_vrml("tri_check", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
/* Gamut center point marker */
wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
@@ -12166,6 +12790,8 @@ double xv[MXRO] /* Intersection point */
/* Main summary plot at each thinning round and at end. */
/* Show vertex surface & optional added or deleted vertexes, */
/* + optional bxcells. */
+#define VV vv /* Actual surface values */
+//#define VV vl /* Logf mapped surface values */
static void plot_vtx_surface(
rspl *s,
int dovtxlabels, /* Show vertex index numbers */
@@ -12193,9 +12819,9 @@ assdire *edgdir /* Edge lookup for vertex */
bxcell *vbx;
if (dopres)
- wrl = new_vrml("final_surface", 0, vrml_lab);
+ wrl = new_vrml("last_surface", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
else
- wrl = new_vrml("thinned_surface", 0, vrml_lab);
+ wrl = new_vrml("thinned_surface", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
if (dovtxlabels) {
@@ -12209,7 +12835,7 @@ assdire *edgdir /* Edge lookup for vertex */
|| (dopres && vx->pres)
|| (dooil && vx->status == vtx_oil)) {
sprintf(index, "%d",vx->ix);
- wrl->add_text(wrl, index, vx->v, cyan, 0.3);
+ wrl->add_text(wrl, index, vx->VV, cyan, 0.3);
}
}
}
@@ -12218,20 +12844,25 @@ assdire *edgdir /* Edge lookup for vertex */
/* Go through the vertex hash to set every vertex value */
for (i = 0; i < vc->hash_size; i++) {
for (vx = vc->hash[i]; vx != NULL; vx = vx->hlink) {
+ double *col = NULL;
if (vx->status != vtx_norm && vx->addvtx)
error ("Found vertex that is both deleted and cause of added bxcell");
if (doadded && vx->addvtx) /* Cause of added bxcell */
- vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, green);
+ col = green;
else if (dopres && vx->pres) /* Preserved vertex */
- vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, yellow);
+ col = yellow;
else if (dodeleted && (vx->status == vtx_sha || vx->status == vtx_del))
- vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, red);
+ col = red;
else if (dooil && vx->status == vtx_oil)
- vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, blue);
+ col = blue;
else if (vx->status == vtx_norm)
- vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->v, white);
+ col = white;
+
+ if (col != NULL) {
+ vx->vrmlix = wrl->add_col_vertex(wrl, 0, vx->VV, col);
+ }
}
}
@@ -12304,26 +12935,26 @@ assdire *edgdir /* Edge lookup for vertex */
col = white;
for (f = 0; f < fdi; f++)
- vv[f] = vx->v[f] + off;
+ vv[f] = vx->VV[f] + off;
lix[0] = wrl->add_vertex(wrl, 2, vv);
for (f = 0; f < fdi; f++)
- vv[f] = vx->v[f] - off;
+ vv[f] = vx->VV[f] - off;
lix[1] = wrl->add_vertex(wrl, 2, vv);
wrl->add_col_line(wrl, 2, lix, col);
for (f = 0; f < fdi; f++)
- vv[f] = vx->v[f] + ((f & 1) ? off : -off);
+ vv[f] = vx->VV[f] + ((f & 1) ? off : -off);
lix[0] = wrl->add_vertex(wrl, 2, vv);
for (f = 0; f < fdi; f++)
- vv[f] = vx->v[f] - ((f & 1) ? off : -off);
+ vv[f] = vx->VV[f] - ((f & 1) ? off : -off);
lix[1] = wrl->add_vertex(wrl, 2, vv);
wrl->add_col_line(wrl, 2, lix, col);
for (f = 0; f < fdi; f++)
- vv[f] = vx->v[f] + ((f & 2) ? off : -off);
+ vv[f] = vx->VV[f] + ((f & 2) ? off : -off);
lix[0] = wrl->add_vertex(wrl, 2, vv);
for (f = 0; f < fdi; f++)
- vv[f] = vx->v[f] - ((f & 2) ? off : -off);
+ vv[f] = vx->VV[f] - ((f & 2) ? off : -off);
lix[1] = wrl->add_vertex(wrl, 2, vv);
wrl->add_col_line(wrl, 2, lix, col);
}
@@ -12396,7 +13027,7 @@ int bxix /* Index of bx cell causing touches */
double white[3] = { 1.0, 1.0, 1.0 };
double red[3] = { 0.8, 0.1, 0.1 };
- wrl = new_vrml("add_touch_bxcells", 0, vrml_lab);
+ wrl = new_vrml("add_touch_bxcells", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
/* Gamut center point marker */
wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
@@ -12461,7 +13092,7 @@ int dowait /* Wait for a return key */
double grey[3] = { 0.5, 0.5, 0.5 };
double white[3] = { 1.0, 1.0, 1.0 };
- wrl = new_vrml("thinned_fwcells", 0, vrml_lab);
+ wrl = new_vrml("thinned_fwcells", 0, s->rev.probxyz ? vrml_xyz : vrml_lab);
wrl->add_marker(wrl, s->rev.ocent, NULL, 1.0);
if (dofclabels) {
diff --git a/rspl/rev.h b/rspl/rev.h
index 6fc79cb..80b9fa9 100644..100755
--- a/rspl/rev.h
+++ b/rspl/rev.h
@@ -62,10 +62,10 @@
* Reversing the Parameter -> Baricentric equations gives the
* following interpolation equation using Parameter coordinates:
*
- * VV0 - VV1 * P0
- * + VV1 - VV2 * P1
- * + VV2 - VV3 * P2
- * + VV3
+ * (VV0 - VV1) * P0
+ * + (VV1 - VV2) * P1
+ * + (VV2 - VV3) * P2
+ * + VV3
*
* Note that withing the simplex, 0 <= P0 && P0 <= P1 && P1 <= P2 && P2 <= 1
*
@@ -270,6 +270,7 @@ struct _bxcell{
bxstat status; /* State of sl list */
int *sl; /* fwd vertex seed list for surface bxcells */
+ /* or cell list after conversion to cells */
int *dl; /* deleted fwd vertex list for this bxcell */
int *scell; /* If non-NULL, this is a (non-surface) */
@@ -441,6 +442,8 @@ struct _rev_struct {
/* All other sections depend on this. */
int fastsetup; /* Flag - NZ if fast setup at cost of slow throughput */
+ int probxyz; /* Flag - NZ - guess if XYZ for VRML diagnostics use */
+
int lchweighted; /* Non-zero if nearest search is LCh weighted */
double lchw[MXRO]; /* LCh nearest weighting */
double lchw_sq[MXRO]; /* LCh nearest weighting squared */
@@ -465,7 +468,7 @@ struct _rev_struct {
int rev_valid; /* nz if this information in rev[] and nnrev[] is valid */
int **rev; /* Exact reverse lookup starting list. */
int **nnrev; /* Nearest reverse lookup starting list. */
- /* These lists are of fwd grid indexes. */
+ /* These lists are of fwd grid base indexes. */
/* [0] is allocation length */
/* [1] is the next free entry index (length + 3, not counting -1) */
/* [2] is index into share lists, -1 if not shared. */
@@ -474,6 +477,10 @@ struct _rev_struct {
double ocent[MXRO]; /* rev cell gamut "center" point for thinning and shadow testing. */
+ int surflin_en; /* Flag set when suflin is enabled */
+ struct _rspl *surflin; /* gamut surface linearization transform used by logcomp() */
+ double linoff[MXRO]; /* ocent offset after surflin mapping */
+
bxcell *surflist; /* Linked list of rev[] bwd cells that contain gamut surface fwd cells. */
/* Used to speed up fill_nncell() when rev.fastsetup is set, else NULL */
int surf_hash_size; /* Current size of bxcell hash list */
diff --git a/rspl/revbench.c b/rspl/revbench.c
index f61ff2e..f61ff2e 100644..100755
--- a/rspl/revbench.c
+++ b/rspl/revbench.c
diff --git a/rspl/rspl.c b/rspl/rspl.c
index db63a47..db63a47 100644..100755
--- a/rspl/rspl.c
+++ b/rspl/rspl.c
diff --git a/rspl/rspl.h b/rspl/rspl.h
index 7d096f9..dc4e14d 100644..100755
--- a/rspl/rspl.h
+++ b/rspl/rspl.h
@@ -24,13 +24,15 @@
/** General Limits **/
-#define MXDI 10 /* Maximum input dimensionality */
-#define MXDO 10 /* Maximum output dimensionality (Is not fully tested!!!) */
-#define LOG2MXDI 4 /* log2 MXDI */
-#define DEF2MXDI 16 /* Default allocation size for 2^di (di=4) */
-#define POW2MXDI 1024 /* 2 ^ MXDI */
-#define DEF3MXDI 81 /* Default allocation size for 3^di (di=4) */
-#define POW3MXDI 59049 /* 3 ^ MXDI */
+#define MXDI 10 /* Maximum input dimensionality */
+#define MXDO 10 /* Maximum output dimensionality (Is not fully tested!!!) */
+#define LOG2MXDI 4 /* log2 MXDI */
+#define DEF2MXDI 16 /* Default allocation size for 2^di (di=4) */
+#define POW2MXDI 1024 /* 2 ^ MXDI */
+#define DEF3MXDI 81 /* Default allocation size for 3^di (di=4) */
+#define POW3MXDI 59049 /* 3 ^ MXDI */
+#define HACOMPS ((POW3MXDI + 2 * MXDI + 1)/2) /* Maximum number of array components */
+ /* 295255 - not used ? */
#if MXDI > MXDO /* Maximum of either DI or DO */
# define MXDIDO MXDI
@@ -38,14 +40,13 @@
# define MXDIDO MXDO
#endif
-/* RESTRICTED SIZE Limits, used for reverse, spline and scattered interpolation */
+/* RESTRICTED SIZE Limits, used for reverse and spline */
#define MXRI 4 /* Maximum input dimensionality */
#define MXRO 10 /* Maximum output dimensionality (Is not fully tested!!!) */
#define LOG2MXRI 2 /* log2 MXRI */
#define POW2MXRI 16 /* 2 ^ MXRI */
#define POW3MXRI 81 /* 3 ^ MXRI */
-#define HACOMPS ((POW3MXRI + 2 * MXRI + 1)/2) /* Maximum number of array components */
#define POW2MXRO 1024 /* 2 ^ MXRO */
#if MXRI > MXRO /* Maximum of either RI or RO */
@@ -93,9 +94,9 @@ typedef struct {
/* Scattered data Per data point data (internal) */
struct _rpnts {
- double p[MXRI]; /* Data position [di] */
- double v[MXRO]; /* Data value [fdi] */
- double k[MXRO]; /* Weight factor (nominally 1.0, less for lower confidence data point) */
+ double p[MXDI]; /* Data position [di] */
+ double v[MXDO]; /* Data value [fdi] */
+ double k[MXDO]; /* Weight factor (nominally 1.0, less for lower confidence data point) */
// double fe; /* Fit error in output pass (ausm) */
}; typedef struct _rpnts rpnts;
@@ -114,7 +115,7 @@ typedef struct {
typedef struct {
int niters; /* Number of multigrid itterations needed */
int **ires; /* Resolution for each itteration and dimension */
- void **mgtmps[MXRO]; /* Store pointers to re-usable mgtmp when incremental */
+ void **mgtmps[MXDO]; /* Store pointers to re-usable mgtmp when incremental */
/* (These don't seem to be used anymore. was incremental removed ?) */
} it_info;
@@ -260,7 +261,7 @@ struct _rspl {
#define RSPL_VERBOSE 0x8000 /* Turn on print progress messages */
#define RSPL_NOVERBOSE 0x4000 /* Turn off print progress messages */
- /* Initialise from scattered data. RESTRICTED SIZE */
+ /* Initialise from scattered data. */
/* Return non-zero if result is non-monotonic */
int
(*fit_rspl)(
@@ -282,7 +283,7 @@ struct _rspl {
/* gres[] entries per dimension. Used to scale smoothness criteria */
);
- /* Initialise from scattered data, with per point weighting. RESTRICTED SIZE */
+ /* Initialise from scattered data, with per point weighting. */
/* Return non-zero if result is non-monotonic */
int
(*fit_rspl_w)(
@@ -305,7 +306,7 @@ struct _rspl {
);
/* Initialise from scattered data, with per point individual out weighting. */
- /* RESTRICTED SIZE Return non-zero if result is non-monotonic */
+ /* Return non-zero if result is non-monotonic */
int
(*fit_rspl_ww)(
struct _rspl *s, /* this */
@@ -327,7 +328,6 @@ struct _rspl {
);
/* Initialise from scattered data, with weak default function. */
- /* RESTRICTED SIZE */
/* Return non-zero if result is non-monotonic */
int
(*fit_rspl_df)(
@@ -353,7 +353,6 @@ struct _rspl {
);
/* Initialise from scattered data, with per point weighting and weak default function. */
- /* RESTRICTED SIZE */
/* Return non-zero if result is non-monotonic */
int
(*fit_rspl_w_df)(
diff --git a/rspl/rspl1.c b/rspl/rspl1.c
index 3996269..3996269 100644..100755
--- a/rspl/rspl1.c
+++ b/rspl/rspl1.c
diff --git a/rspl/rspl1.h b/rspl/rspl1.h
index 63be428..63be428 100644..100755
--- a/rspl/rspl1.h
+++ b/rspl/rspl1.h
diff --git a/rspl/rspl_imp.h b/rspl/rspl_imp.h
index 0cd8805..0cd8805 100644..100755
--- a/rspl/rspl_imp.h
+++ b/rspl/rspl_imp.h
diff --git a/rspl/scat.c b/rspl/scat.c
index f94e074..87835ea 100644..100755
--- a/rspl/scat.c
+++ b/rspl/scat.c
@@ -221,20 +221,24 @@ struct _mgtmp {
int bres, brix; /* Biggest resolution and its index */
double mres; /* Geometric mean res[] */
int no; /* Total number of points in grid = res ^ di */
- ratai l,h,w; /* Grid low, high, grid cell width */
+ datai l,h,w; /* Grid low, high, grid cell width */
double *ipos[MXDI]; /* Optional relative grid cell position for each input dim cell */
/* Grid array offset lookups */
- int ci[MXRI]; /* Grid coordinate increments for each dimension */
- int hi[POW2MXRI]; /* Combination offset for sequence through cube. */
+ int ci[MXDI]; /* Grid coordinate increments for each dimension */
+ int hi[POW2MXDI]; /* Combination offset for sequence through cube. */
} g;
- /* Data point grid dependent information */
+ /* Data point grid dependent information - allocate for size of 2^di */
struct mgdat {
int b; /* Index for associated base grid point, in grid points */
- double w[POW2MXRI]; /* Weight for surrounding gridpoints [2^di] */
+ double w[0]; /* Weight for surrounding gridpoints [2^di] */
} *d;
+ int sizeof_mgdat; /* sizeof(mgdat) + sizeof(double) * (1<<di) */
+
+ /* Address of mgdat at index N from given mgtmp *M */
+#define MGDAT_N(M,N) ((struct mgdat *)(((char *)&(M)->d[0]) + (N) * (M)->sizeof_mgdat))
/* Equation Solution related (Grid point solutions) */
struct {
@@ -245,8 +249,7 @@ struct _mgtmp {
/* and +/-1 offset in all dimensions, but only the +ve offset */
/* half of the sparse matrix is stored, due to equations */
/* being symetrical. */
-#define XCOLPMAX (HACOMPS+8)
- int xcol[XCOLPMAX]; /* A array column translation from packed to sparse index */
+ int *xcol; /* A array column translation from packed to sparse index */
int *ixcol; /* A array column translation from sparse to packed index */
double *b; /* b vector for RHS of simultabeous equation b[g.no] */
double normb; /* normal of b vector */
@@ -408,11 +411,11 @@ fit_rspl_imp(
void *d, /* Array holding position and function values of data points */
int dtp, /* Flag indicating data type, 0 = (co *), 1 = (cow *), 2 = (coww *) */
int dno, /* Number of data points */
- ratai glow, /* Grid low scale - will be expanded to enclose data, NULL = default 0.0 */
- ratai ghigh, /* Grid high scale - will be expanded to enclose data, NULL = default 1.0 */
+ 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[MXDI], /* Spline grid resolution */
- ratao vlow, /* Data value low normalize, NULL = default 0.0 */
- ratao vhigh, /* Data value high normalize - NULL = default 1.0 */
+ datao vlow, /* Data value low normalize, NULL = default 0.0 */
+ datao vhigh, /* Data value high normalize - NULL = default 1.0 */
double smooth, /* Smoothing factor, 0.0 = default 1.0 */
/* (if -ve, overides optimised smoothing, and sets raw smoothing */
/* typically between 1e-7 .. 1e-1) */
@@ -430,12 +433,11 @@ fit_rspl_imp(
int i, e, f;
#ifdef NEVER
-printf("~1 rspl: gres = %d %d %d %d, smooth = %f, avgdev = %f %f %f\n",
-gres[0], gres[1], gres[2], gres[3], smooth, avgdev[0], avgdev[1], avgdev[2]);
-printf("~1 rspl: glow = %f %f %f %f ghigh = %f %f %f %f\n",
-glow[0], glow[1], glow[2], glow[3], ghigh[0], ghigh[1], ghigh[2], ghigh[3]);
-printf("~1 rspl: vlow = %f %f %f vhigh = %f %f %f\n",
-vlow[0], vlow[1], vlow[2], vhigh[0], vhigh[1], vhigh[2]);
+printf("~1 rspl: di %d, fdi %d\n",fdi);
+printf("~1 rspl: gres = %s, smooth = %f, avgdev = %s\n",
+debPiv(di, gres),smooth, debPdv(fdi,avgdev));
+printf("~1 rspl: glow = %s ghigh = %s\n", debPdv(di,glow), debPdv(di,ghigh));
+printf("~1 rspl: vlow = %s vhigh = %s\n",debPdv(fdi,vlow), debPdv(fdi,vhigh));
printf("~1 rspl: flags = 0x%x\n",flags);
#endif
@@ -444,9 +446,9 @@ printf("~1 rspl: flags = 0x%x\n",flags);
#endif
/* This is a restricted size function */
- if (di > MXRI)
+ if (di > MXDI)
error("rspl: fit can't handle di = %d",di);
- if (fdi > MXRO)
+ if (fdi > MXDO)
error("rspl: fit can't handle fdi = %d",fdi);
/* set debug level */
@@ -528,8 +530,9 @@ printf("~1 rspl: flags = 0x%x\n",flags);
for (i = 0; i < dno; i++) {
for (e = 0; e < s->di; e++) {
- if (dp[i].p[e] > s->g.h[e])
+ if (dp[i].p[e] > s->g.h[e]) {
s->g.h[e] = dp[i].p[e];
+ }
if (dp[i].p[e] < s->g.l[e])
s->g.l[e] = dp[i].p[e];
}
@@ -882,11 +885,11 @@ fit_rspl(
int flags, /* Combination of flags */
co *d, /* Array holding position and function values of data points */
int dno, /* Number of data points */
- ratai glow, /* Grid low scale - will be expanded to enclose data, NULL = default 0.0 */
- ratai ghigh, /* Grid high scale - will be expanded to enclose data, NULL = default 1.0 */
+ 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[MXDI], /* Spline grid resolution */
- ratao vlow, /* Data value low normalize, NULL = default 0.0 */
- ratao vhigh, /* Data value high normalize - NULL = default 1.0 */
+ 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 */
double avgdev[MXDO],
/* Average Deviation of function values as proportion of function range. */
@@ -906,11 +909,11 @@ fit_rspl_w(
int flags, /* Combination of flags */
cow *d, /* Array holding position, function and weight values of data points */
int dno, /* Number of data points */
- ratai glow, /* Grid low scale - will be expanded to enclose data, NULL = default 0.0 */
- ratai ghigh, /* Grid high scale - will be expanded to enclose data, NULL = default 1.0 */
+ 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[MXDI], /* Spline grid resolution */
- ratao vlow, /* Data value low normalize, NULL = default 0.0 */
- ratao vhigh, /* Data value high normalize - NULL = default 1.0 */
+ 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 */
double avgdev[MXDO],
/* Average Deviation of function values as proportion of function range. */
@@ -930,11 +933,11 @@ fit_rspl_ww(
int flags, /* Combination of flags */
coww *d, /* Array holding position, function and weight values of data points */
int dno, /* Number of data points */
- ratai glow, /* Grid low scale - will be expanded to enclose data, NULL = default 0.0 */
- ratai ghigh, /* Grid high scale - will be expanded to enclose data, NULL = default 1.0 */
+ 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[MXDI], /* Spline grid resolution */
- ratao vlow, /* Data value low normalize, NULL = default 0.0 */
- ratao vhigh, /* Data value high normalize - NULL = default 1.0 */
+ 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 */
double avgdev[MXDO],
/* Average Deviation of function values as proportion of function range. */
@@ -1530,15 +1533,17 @@ static mgtmp *new_mgtmp(
#endif
/* Allocate space for auiliary data point related info */
- if ((m->d = (struct mgdat *) calloc(dno, sizeof(struct mgdat))) == NULL)
+ m->sizeof_mgdat = sizeof(struct mgdat) + sizeof(double) * (1 << di);
+ if ((m->d = (struct mgdat *) calloc(dno, m->sizeof_mgdat)) == NULL)
error("rspl: malloc failed - mgtmp");
/* fill in the aux data point info */
/* (We're assuming N-linear interpolation here. */
/* Perhaps we should try simplex too ?) */
for (n = 0; n < dno; n++) {
- double we[MXRI]; /* 1.0 - Weight in each dimension */
+ double we[MXDI]; /* 1.0 - Weight in each dimension */
int ix = 0; /* Index to base corner of surrounding cube in grid points */
+ struct mgdat *dd;
/* Figure out which grid cell the point falls into */
for (e = 0; e < di; e++) {
@@ -1557,21 +1562,22 @@ static mgtmp *new_mgtmp(
ix += mi * m->g.ci[e]; /* Add Index offset for grid cube base in dimen */
we[e] = t - (double)mi; /* 1.0 - weight */
}
- m->d[n].b = ix;
+ dd = MGDAT_N(m, n);
+ dd->b = ix;
/* Compute corner weights needed for interpolation */
- m->d[n].w[0] = 1.0;
+ dd->w[0] = 1.0;
for (e = 0, g = 1; e < di; g *= 2, e++) {
for (i = 0; i < g; i++) {
- m->d[n].w[g+i] = m->d[n].w[i] * we[e];
- m->d[n].w[i] *= (1.0 - we[e]);
+ dd->w[g+i] = dd->w[i] * we[e];
+ dd->w[i] *= (1.0 - we[e]);
}
}
#ifdef DEBUG
printf("Data point %d weighting factors = \n",n);
for (e = 0; e < (1 << di); e++) {
- printf("%d: %f\n",e,m->d[n].w[e]);
+ printf("%d: %f\n",e,dd->w[e]);
}
#endif /* DEBUG */
@@ -1588,6 +1594,7 @@ static mgtmp *new_mgtmp(
/* Set the solution matricies to unalocated */
m->q.A = NULL;
+ m->q.xcol = NULL;
m->q.ixcol = NULL;
m->q.b = NULL;
m->q.x = NULL;
@@ -1620,6 +1627,7 @@ static void free_mgtmp(mgtmp *m) {
}
free_dvector(m->q.x,0,gno-1);
free_dvector(m->q.b,0,gno-1);
+ free((void *)m->q.xcol);
free((void *)m->q.ixcol);
free_dmatrix(m->q.A,0,gno-1,0,m->q.acols-1);
free((void *)m->d);
@@ -1770,8 +1778,24 @@ mgtmp *sm /* Optional smoothing map for ausm mode */
DCOUNT(gc, MXDIDO, di, -2, -2, 3); /* Step through +/- 2 cube offset */
#endif
int ix; /* Grid point offset in grid points */
+ int xcolpmax;
acols = 0;
+ /* Allocate xcol[] */
+ if (xcol == NULL) {
+ /* Comput used dimensions XCOLPMAX, */
+ /* = number of array components, */
+ /* = HACOMPS + 8 */
+ xcolpmax = 1;
+ for (k = 0; k < di; k++)
+ xcolpmax *= 3;
+ xcolpmax += 2 * di + 1;
+ xcolpmax /= 2;
+ xcolpmax += 8;
+ if ((xcol = (int *) malloc(xcolpmax * sizeof(int))) == NULL)
+ error("rspl malloc failed - xcol");
+ }
+
DC_INIT(gc);
while (!DC_DONE(gc)) {
@@ -1796,7 +1820,7 @@ mgtmp *sm /* Optional smoothing map for ausm mode */
for (ix = 0, k = 0; k < di; k++)
ix += gc[k] * gci[k]; /* Multi-dimension grid offset */
if (ix >= 0) {
- if (acols >= XCOLPMAX)
+ if (acols >= xcolpmax)
error("rspl internal: exceeded xcol bounds");
xcol[acols++] = ix; /* We only store half, due to symetry */
}
@@ -1849,6 +1873,7 @@ mgtmp *sm /* Optional smoothing map for ausm mode */
m->q.b = b;
m->q.x = x;
m->q.acols = acols;
+ m->q.xcol = xcol;
m->q.ixcol = ixcol;
} else { /* re-initializing, zero matrices */
@@ -2264,7 +2289,8 @@ mgtmp *sm /* Optional smoothing map for ausm mode */
/* Accumulate data point dependent factors */
for (n = 0; n < dno; n++) { /* Go through all the data points */
int j,k;
- int bp = m->d[n].b; /* index to base grid point in grid points */
+ struct mgdat *dd = MGDAT_N(m, n);
+ int bp = dd->b; /* index to base grid point in grid points */
/* For each point in the cube as the base grid point, */
/* add in the appropriate weighting for its weighted neighbors. */
@@ -2274,7 +2300,7 @@ mgtmp *sm /* Optional smoothing map for ausm mode */
ai = bp + m->g.hi[j]; /* A matrix index */
- w = m->d[n].w[j]; /* Base point grid weight */
+ w = dd->w[j]; /* Base point grid weight */
d = 2.0 * s->d.a[n].k[f] * w; /* (2.0, w are derivtv factors, k data pnt wgt) */
tt = d * s->d.a[n].v[f]; /* Change in data component */
@@ -2287,7 +2313,7 @@ mgtmp *sm /* Optional smoothing map for ausm mode */
for (k = j+1; k < (1 << di); k++) { /* Binary sequence */
int ii;
ii = ixcol[m->g.hi[k] - m->g.hi[j]]; /* A matrix column index */
- A[ai][ii] += d * m->d[n].w[k]; /* dui component due to ui+1 */
+ A[ai][ii] += d * dd->w[k]; /* dui component due to ui+1 */
}
}
}
@@ -2337,14 +2363,15 @@ static void comp_fit_errors(
/* Compute error for each data point */
for (n = 0; n < dno; n++) {
int j;
- int bp = m->d[n].b; /* index to base grid point in grid points */
+ struct mgdat *dd = MGDAT_N(m, n);
+ int bp = dd->b; /* index to base grid point in grid points */
double val; /* Current interpolated value */
double err;
double gain = 1.0;
/* Compute the interpolated grid value for this data point */
for (val = 0.0, j = 0; j < (1 << di); j++) /* Binary sequence */
- val += m->d[n].w[j] * x[bp + m->g.hi[j]];
+ val += dd->w[j] * x[bp + m->g.hi[j]];
err = s->d.a[n].v[f] - val;
err *= 0.8;
@@ -2367,8 +2394,8 @@ static double mgtmp_interp(
rspl *s = m->s;
int di = s->di;
int e, n;
- double we[MXRI]; /* 1.0 - Weight in each dimension */
- double gw[POW2MXRI]; /* weight for each grid cube corner */
+ double we[MXDI]; /* 1.0 - Weight in each dimension */
+ double gw[POW2MXDI]; /* weight for each grid cube corner */
double *gp; /* Pointer to x2[] grid cube base */
double val;
@@ -2428,7 +2455,7 @@ static void init_soln(
/* For all output grid points */
EC_INIT(gc);
for (n = 0; n < gno; n++) {
- double p[MXRI]; /* Grid relative location */
+ double p[MXDI]; /* Grid relative location */
for (e = 0; e < di; e++)
p[e] = (double)gc[e]/(m1->g.res[e] - 1.0);
@@ -2645,7 +2672,7 @@ one_itter1(
/* For each dimension */
for (d = 0; d < di; d++) {
int ld = d == 0 ? 1 : 0; /* lowest dim */
- int sof, gc[MXRI];
+ int sof, gc[MXDI];
//printf("~1 doing one_itter1 for dim %d\n",d);
for (e = 0; e < di; e++)
@@ -2713,7 +2740,7 @@ one_itter2(
double ovsh /* Overshoot to use, 1.0 for none */
) {
int e,i,k;
- int gc[MXRI];
+ int gc[MXDI];
for (i = e = 0; e < di; e++)
gc[e] = 0; /* init coords */
diff --git a/rspl/scat2.c b/rspl/scat2.c
index 7579366..7579366 100644..100755
--- a/rspl/scat2.c
+++ b/rspl/scat2.c
diff --git a/rspl/sm1.c b/rspl/sm1.c
index 5bcf26e..5bcf26e 100644..100755
--- a/rspl/sm1.c
+++ b/rspl/sm1.c
diff --git a/rspl/sm2.c b/rspl/sm2.c
index 84bc8cf..84bc8cf 100644..100755
--- a/rspl/sm2.c
+++ b/rspl/sm2.c
diff --git a/rspl/sm3.c b/rspl/sm3.c
index 9307787..9307787 100644..100755
--- a/rspl/sm3.c
+++ b/rspl/sm3.c
diff --git a/rspl/smtmpp.c b/rspl/smtmpp.c
index 1264d4b..1264d4b 100644..100755
--- a/rspl/smtmpp.c
+++ b/rspl/smtmpp.c
diff --git a/rspl/smtnd.c b/rspl/smtnd.c
index f105429..355d20b 100644..100755
--- a/rspl/smtnd.c
+++ b/rspl/smtnd.c
@@ -1031,7 +1031,7 @@ static void do_test(
if (verb) printf("Measured noise average deviation = %f%%\n",tnoise * 100.0);
/* Fit to scattered data */
- if (verb) printf("Fitting the scattered data, smooth = %f, avgdev = %f\n",smooth,avgdev != NULL ? avgdev[0] : 0.0);
+ if (verb) printf("Fitting the scattered data, smooth = %f, avgdev = %f\n",smooth, avgdev[0]);
avgdev[0] = 0.25 * noise;
rss->fit_rspl(rss,
flags, /* Non-mon and clip flags */
diff --git a/rspl/spline.c b/rspl/spline.c
index dabeb0f..dabeb0f 100644..100755
--- a/rspl/spline.c
+++ b/rspl/spline.c
diff --git a/rspl/stest.c b/rspl/stest.c
index d1cf2da..d1cf2da 100644..100755
--- a/rspl/stest.c
+++ b/rspl/stest.c
diff --git a/rspl/t2d.c b/rspl/t2d.c
index e3f304c..e3f304c 100644..100755
--- a/rspl/t2d.c
+++ b/rspl/t2d.c
diff --git a/rspl/t2ddf.c b/rspl/t2ddf.c
index 115d5f5..115d5f5 100644..100755
--- a/rspl/t2ddf.c
+++ b/rspl/t2ddf.c
diff --git a/rspl/t3d.c b/rspl/t3d.c
index c526669..c526669 100644..100755
--- a/rspl/t3d.c
+++ b/rspl/t3d.c
diff --git a/rspl/t3ddf.c b/rspl/t3ddf.c
index 4a43f0d..4a43f0d 100644..100755
--- a/rspl/t3ddf.c
+++ b/rspl/t3ddf.c
diff --git a/rspl/tnd.c b/rspl/tnd.c
index 89718d9..89718d9 100644..100755
--- a/rspl/tnd.c
+++ b/rspl/tnd.c
diff --git a/rspl/trnd.c b/rspl/trnd.c
index 2676182..2676182 100644..100755
--- a/rspl/trnd.c
+++ b/rspl/trnd.c
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-27 05:26:22 +0000