406 lines
9.0 KiB
C
406 lines
9.0 KiB
C
|
|
/**********************************************************/
|
|
/* 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>
|
|
#include "aconfig.h"
|
|
#include "numlib.h"
|
|
#include "plot.h"
|
|
#include "ui.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";
|
|
|
|
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");
|
|
}
|
|
|