summaryrefslogtreecommitdiff
path: root/xicc/cvtest.c
diff options
context:
space:
mode:
Diffstat (limited to 'xicc/cvtest.c')
-rw-r--r--xicc/cvtest.c411
1 files changed, 411 insertions, 0 deletions
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 <stdio.h>
+#include <stdlib.h>
+#include <fcntl.h>
+#include <math.h>
+#if defined(__IBMC__) && defined(_M_IX86)
+#include <float.h>
+#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");
+}
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-10-31 23:05:02 +0000