diff options
Diffstat (limited to 'rspl')
| -rw-r--r-- | rspl/Jamfile | 8 | ||||
| -rw-r--r-- | rspl/afiles | 2 | ||||
| -rw-r--r-- | rspl/c1.c | 3 | ||||
| -rw-r--r-- | rspl/cw1.c | 308 | ||||
| -rw-r--r-- | rspl/cw3.c | 325 | ||||
| -rw-r--r-- | rspl/rev.c | 34 | ||||
| -rw-r--r-- | rspl/rspl.c | 2 | ||||
| -rw-r--r-- | rspl/scat.c | 270 |
8 files changed, 752 insertions, 200 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/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 @@ -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 */ @@ -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++) { |