From 22f703cab05b7cd368f4de9e03991b7664dc5022 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Mon, 1 Sep 2014 13:56:46 +0200 Subject: Initial import of argyll version 1.5.1-8 --- xicc/cvtest.c | 411 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 411 insertions(+) create mode 100644 xicc/cvtest.c (limited to 'xicc/cvtest.c') diff --git a/xicc/cvtest.c b/xicc/cvtest.c new file mode 100644 index 0000000..6d24572 --- /dev/null +++ b/xicc/cvtest.c @@ -0,0 +1,411 @@ + +/**********************************************************/ +/* Investigate GEMS transfer curve function approximation */ +/**********************************************************/ + +/* Standard test + Random testing */ + +/* Author: Graeme Gill + * Date: 2003/12/1 + * + * Copyright 2003 Graeme W. Gill + * Parts derived from rspl/c1.c + * + * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :- + * see the License.txt file for licencing details. + */ + +#undef DIAG +#undef TEST_SYM /* Test forcing center to zero (a*, b* constraint) */ + +#include +#include +#include +#include +#if defined(__IBMC__) && defined(_M_IX86) +#include +#endif +#include "numlib.h" +#include "plot.h" + +#define MAX_PARM 40 /* Make > SHAPE_ORDS */ + +#define POWTOL 1e-5 +#define MAXITS 10000 + +/* SHAPE_BASE doesn't seem extremely critical. It is centered in +/- 1 magnitude */ +/* 10 x more filters out noise reasonably heaviliy, 10 x less gives noticable */ +/* overshoot Range 0.00001 .. 0.001 */ +/* SHAPE_HBASE is more critical. */ +/* Range 0.00005 .. 0.001 */ +//#define SHAPE_BASE 0.00001 /* 0 & 1 harmonic weight */ +//#define SHAPE_HBASE 0.0002 /* 2nd and higher additional weight */ +//#define SHAPE_ORDS 20 +#define SHAPE_BASE 0.00001 /* 0 & 1 harmonic weight */ +#define SHAPE_HBASE 0.0001 /* 2nd and higher additional weight */ +#define SHAPE_ORDS 20 + +/* Interface coordinate value */ +typedef struct { + double p; /* coordinate position */ + double v; /* function values */ +} co; + +double lin(double x, double xa[], double ya[], int n); +static double tfunc(double *v, int luord, double vv); +void fit(double *params, int np, co *test_points, int pnts); +void usage(void); + +#define TRIALS 40 /* Number of random trials */ +#define SKIP 0 /* Number of random trials to skip */ + +#define ABS_MAX_PNTS 100 + +#define MIN_PNTS 2 +#define MAX_PNTS 20 + +#define MIN_RES 20 +#define MAX_RES 500 + +double xa[ABS_MAX_PNTS]; +double ya[ABS_MAX_PNTS]; + +#define XRES 100 + +#define TSETS 4 +#define PNTS 11 +#define GRES 100 +int t1p[TSETS] = { + 4, + 11, + 11, + 11 +}; + +double t1xa[TSETS][PNTS] = { + { 0.0, 0.2, 0.8, 1.0 }, + { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 }, + { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 }, + { 0.0, 0.25, 0.30, 0.35, 0.40, 0.44, 0.48, 0.51, 0.64, 0.75, 1.0 } +}; + +double t1ya[TSETS][PNTS] = { + { 0.0, 0.5, 0.6, 1.0 }, + { 0.0, 0.10, 0.22, 0.35, 0.52, 0.65, 0.78, 0.91, 1.0, 0.9, 0.85 }, + { 0.0, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.8, 1.0, 1.0, 1.0 }, + { 0.0, 0.35, 0.4, 0.41, 0.42, 0.46, 0.5, 0.575, 0.48, 0.75, 1.0 } +}; + + +co test_points[ABS_MAX_PNTS]; + +double lin(double x, double xa[], double ya[], int n); + +int main(void) { + int i, n; + double x; + double xx[XRES]; + double y1[XRES]; + double y2[XRES]; + int np = SHAPE_ORDS; /* Number of harmonics */ + double params[MAX_PARM]; /* Harmonic parameters */ + + error_program = "Curve1"; + +#if defined(__IBMC__) + _control87(EM_UNDERFLOW, EM_UNDERFLOW); + _control87(EM_OVERFLOW, EM_OVERFLOW); +#endif + + for (n = 0; n < TRIALS; n++) { + double lrand; /* Amount of level randomness */ + int pnts; + + if (n < TSETS) /* Standard versions */ { + pnts = t1p[n]; + for (i = 0; i < pnts; i++) { + xa[i] = t1xa[n][i]; + ya[i] = t1ya[n][i]; + } + } else if (n == TSETS) { /* Exponential function aproximation */ + double ex = 2.2; + pnts = MAX_PNTS; + + printf("Trial %d, no points = %d, exponential %f\n",n,pnts,ex); + + /* Create X values */ + for (i = 0; i < pnts; i++) + xa[i] = i/(pnts-1.0); + + for (i = 0; i < pnts; i++) + ya[i] = pow(xa[i], ex); + + } else if (n == (TSETS+1)) { /* Exponential function aproximation */ + double ex = 1.0/2.2; + pnts = MAX_PNTS; + + printf("Trial %d, no points = %d, exponential %f\n",n,pnts,ex); + + /* Create X values */ + for (i = 0; i < pnts; i++) + xa[i] = i/(pnts-1.0); + + for (i = 0; i < pnts; i++) + ya[i] = pow(xa[i], ex); + + } else { /* Random versions */ + lrand = d_rand(0.0,0.2); /* Amount of level randomness */ + lrand *= lrand; + pnts = i_rand(MIN_PNTS,MAX_PNTS); + + printf("Trial %d, no points = %d, level randomness = %f\n",n,pnts,lrand); + + /* Create X values */ + xa[0] = 0.0; + for (i = 1; i < pnts; i++) + xa[i] = xa[i-1] + d_rand(0.5,1.0); + for (i = 0; i < pnts; i++) /* Divide out */ + xa[i] = (xa[i]/xa[pnts-1]); + + /* Create y values */ + ya[0] = xa[0]; + for (i = 1; i < pnts; i++) + ya[i] = ya[i-1] + d_rand(0.1,1.0) + d_rand(-0.1,0.4) + d_rand(-0.4,0.5); + for (i = 0; i < pnts; i++) { /* Divide out */ + ya[i] = (ya[i]/ya[pnts-1]); + if (ya[i] < 0.0) + ya[i] = 0.0; + else if (ya[i] > 1.0) + ya[i] = 1.0; + } + } + + if (n < SKIP) + continue; + + for (i = 0; i < pnts; i++) { + test_points[i].p = xa[i]; + test_points[i].v = ya[i]; + } + + for (np = 2; np <= SHAPE_ORDS; np++) { + + /* Fit to scattered data */ + fit(params, /* Parameters to return */ + np, /* Number of parameters */ + test_points, /* Test points */ + pnts /* Number of test points */ + ); + + printf("Number params = %d\n",np); + for (i = 0; i < np; i++) { + printf("Param %d = %f\n",i,params[i]); + } + + /* Display the result */ + for (i = 0; i < XRES; i++) { + x = i/(double)(XRES-1); + xx[i] = x; + y1[i] = lin(x,xa,ya,pnts); + y2[i] = tfunc(params, np, x); + if (y2[i] < -0.2) + y2[i] = -0.2; + else if (y2[i] > 1.2) + y2[i] = 1.2; + } + + do_plot(xx,y1,y2,NULL,XRES); + } + + } /* next trial */ + return 0; +} + + +double lin( +double x, +double xa[], +double ya[], +int n) { + int i; + double y; + + if (x < xa[0]) + return ya[0]; + else if (x > xa[n-1]) + return ya[n-1]; + + for (i = 0; i < (n-1); i++) + if (x >=xa[i] && x <= xa[i+1]) + break; + + x = (x - xa[i])/(xa[i+1] - xa[i]); + + y = ya[i] + (ya[i+1] - ya[i]) * x; + + return y; + } + + +/******************************************************************/ +/* Error/debug output routines */ +/******************************************************************/ + +void +usage(void) { + puts("usage: curve"); + exit(1); + } + +/*******************************************************************/ +/* Grapic gems based, monotonic 1D function transfer curve. */ + +/* Per transfer function */ +static double tfunc( +double *v, /* Pointer to first parameter */ +int luord, /* Curve order n..MPP_MXTCORD */ +double vv /* Source of value */ +) { + double g; + int ord; + + if (vv < 0.0) + vv = 0.0; + else if (vv > 1.0) + vv = 1.0; + + /* Process all the shaper orders from high to low. */ + /* [These shapers were inspired by a Graphics Gem idea */ + /* (Gems IV, VI.3, "Fast Alternatives to Perlin's Bias and */ + /* Gain Functions, pp 401). */ + /* They have the nice properties that they are smooth, and */ + /* can't be non-monotonic. The control parameter has been */ + /* altered to have a range from -oo to +oo rather than 0.0 to 1.0 */ + /* so that the search space is less non-linear. ] */ + for (ord = 0; ord < luord; ord++) { + int nsec; /* Number of sections */ + double sec; /* Section */ + + g = v[ord]; /* Parameter */ + + nsec = ord + 1; /* Increase sections for each order */ + + vv *= (double)nsec; + + sec = floor(vv); + if (((int)sec) & 1) + g = -g; /* Alternate action in each section */ + vv -= sec; + if (g >= 0.0) { + vv = vv/(g - g * vv + 1.0); + } else { + vv = (vv - g * vv)/(1.0 - g * vv); + } + vv += sec; + vv /= (double)nsec; + } + + return vv; +} + +/* return the sum of the squares of the current shaper parameters */ +static double shapmag( +double *v, /* Pointer to first parameter */ +int luord +) { + double tt, w, tparam = 0.0; + int f; + + for (f = 0; f < luord; f++) { + tt = v[f]; + tt *= tt; + + /* Weigh to suppress ripples */ + if (f <= 1) + w = SHAPE_BASE; + else { + w = SHAPE_BASE + (f-1) * SHAPE_HBASE; + } + + tparam += w * tt; + } + return tparam; +} + +/* Context for optimising function fit */ +typedef struct { + int luord; /* shaper order */ + co *points; /* List of test points */ + int nodp; /* Number of data points */ +} luopt; + +/* Shaper+Matrix optimisation function handed to powell() */ +static double luoptfunc(void *edata, double *v) { + luopt *p = (luopt *)edata; + double rv = 0.0, smv; + double out; + int i; + + /* For all our data points */ + for (i = 0; i < p->nodp; i++) { + double ev; + + /* Apply our function */ + out = tfunc(v, p->luord, p->points[i].p); + + ev = out - p->points[i].v; + +#ifdef NEVER + ev = fabs(ev); + rv += ev * ev * ev; +#else + rv += ev * ev; +#endif + + } + + /* Normalise error to be an average delta E squared */ + rv /= (double)p->nodp; + + /* Sum with shaper parameters squared, to */ + /* minimise unsconstrained "wiggles" */ + smv = shapmag(v, p->luord); + rv += smv; + +#ifdef TEST_SYM + { + double tt; + tt = tfunc(v, p->luord, 0.5) - 0.5; + tt *= tt; + rv += 200.0 * tt; + } +#endif + + printf("~1 rv = %f (%f)\n",rv,smv); + return rv; +} + +/* Fitting function */ +void fit( +double *params, /* Parameters to return */ +int np, /* Number of parameters */ +co *test_points, /* Test points */ +int pnts /* Number of test points */ +) { + int i; + double sa[MAX_PARM]; /* Search area */ + luopt os; + + os.luord = np; + os.points= test_points; + os.nodp = pnts; + + for (i = 0; i < np; i++) { + sa[i] = 0.5; + params[i] = 0.0; + } + + if (powell(NULL, np, params, sa, POWTOL, MAXITS, luoptfunc, (void *)&os, NULL, NULL) != 0) + error ("Powell failed"); +} + -- cgit v1.2.3