[argyllcms] Randomized GCV for RSPL
- From: Gerhard Fuernkranz <nospam456@xxxxxx>
- To: Graeme Gill <graeme@xxxxxxxxxxxxx>, argyllcms@xxxxxxxxxxxxx
- Date: Sat, 24 Jun 2006 16:50:04 +0200
Hello Graeme,
the attached patch adds the computation of randomized GCV (Generalized
Cross Validation) to the RSPL code. Please take a look, whether you also
find it useful. I've also added a very crude, quick and dirty
optimization routine, which attempts to minimize GCV (and which has
sureley potential for improvements). To get a first impression, please
apply the patch (relative to the current development snapshot), and try
"profile -v -v ..." and "profile -v -O1 ...".
The code is based on
Wahba, G., Johnson, D. R., Gao, F. and Gong, J. "Adaptive tuning
of numerical weather prediction models: Part I: randomized GCV and
related methods in three and four dimensional data assimilation."
TR 920, April 1994.
(http://citeseer.ist.psu.edu/135988.html)
but I used random numbers from a binary distribution [-sigma,sigma]
(each with 50% probablilty) to perturb the data, as proposed by
Hutchinson (http://www.google.de/search?q=hutchinson+trace+estimator),
because my experimental impression was that the latter distribution
gives GCV estimates with a lower variance (I'm not sure, whether this
observation is mathematically founded, or whether it just happend to be
the case in conjunction with the measurement sets I tried).
While a completely analytic calculation of the GCV (and other similar
metrics which attempt to estimate the goodness of the generalization of
a regression model to unkown data) for a generalized linear model would
require the inversion of the influence matrix, which is numerically
rather no longer feasible for a spline model, due to the size of the
problem (e.g. 35937x35937 matrix for a 33-point 3D grid), this
stochastic method allows the estimation of GCV (and also of other
similar metrics) by fitting the regression model for a 2nd time to a
randomly perturbed variant of the training set. The method has also the
nice property that it treats the regression model as a black box, and it
can be even applied approximately to non-linear models, if the magnitude
of the random noise used to perturb the data is chosen carefully. But it
is a stochastic method, and thus the resulting GCV estimate has some
variance, i.e. each different vector of random numbers used to perturb
the traning set will result in a slightly different GCV estimate. The
larger the training set, the lower the variance will become; for small
training sets, the trace estimates resulting from several different
random perturbation vectors can be averaged to lower the variance (but
such an averaging requires to fit the data not just for a 2nd time, but
several times, imposing a higher computational cost). To get more
consist results when comparing the GCV estimates for different models
(e.g. with different smoothness factors), it is helpful to use the same
random vector(s) for the GCV computation for each of the compared models
(as it is also suggested in the paper).
One of my observations with a couple of data sets is furthermore, that
the three output dimensions tend to require a different amount of
smoothness. Depending on the data set being fitted (i.e. depending on
the characteristics of the device being profiled), the L* channel mostly
seems to require more smoothness, typically by a factor of 3..30 higher
than for the a* and b* channels. Particularly for scanner data sets, I
encountered larger optimal smoothing factor values for L*. So it may
probably also make sense, if fit_rspl() and fiends would take not just a
scalar "smooth" argument, but if a vector could be supplied.
Regards,
Gerhard
*** rspl/rspl.h.orig 2006-06-23 23:00:37.000000000 +0200
--- rspl/rspl.h 2006-06-23 20:56:45.000000000 +0200
***************
*** 204,209 ****
--- 204,218 ----
/* Reverse Interpolation support */
rev_struct rev; /* See rev.h */
+ /* GCV support */
+ struct rspl_gcv {
+ double rss[MXDO]; /* residual sum of quares for each
output channel */
+ double df[MXDO]; /* degrees of freedom */
+ double effpar[MXDO]; /* effective number of parameters */
+ double gcv[MXDO]; /* sigma^2 GCV estimate for each output
channel */
+ double gcv_sum; /* overall sigma^2 GCV estimate for all
channels */
+ } gcv;
+
/* Methods */
/* Free ourselves */
***************
*** 215,220 ****
--- 224,231 ----
#define RSPL_SYMDOMAIN 0x0004 /* Maintain symetric smoothness
with nonsym. resolution */
#define RSPL_INCREMENTAL 0x0008 /* Enable adding more data
points */
#define RSPL_FINAL 0x0010 /* Signal to add_rspl() that
this is the last points */
+ #define RSPL_GCV 0x0020 /* Compute randomized GCV
estimate */
+ #define RSPL_OPT_GCV 0x0040 /* Optimize smoothing for
minimal GCV */
#define RSPL_VERBOSE 0x8000 /* Print progress messages */
/* Initialise from scattered data. RESTRICTED SIZE */
*** rspl/scat.c.orig 2006-05-26 06:06:25.000000000 +0200
--- rspl/scat.c 2006-06-23 23:43:33.000000000 +0200
***************
*** 132,142 ****
--- 132,150 ----
#define NME_TARGET 0.5 /* nme error target */
+ #define GCV_MIN_POINTS 250 /* Minimum number of data points used for the
+ randomized trace estimataion. The higher, the
+ lower the variance of the trace estimate. But
+ a high number will cost additional
iterations */
+ /* #define GCV_MIN_POINTS 1000 */
+
#undef DEBUG_PROGRESS /* Print progress of acheiving tollerance target */
#undef NEVER
#define ALWAYS
+ #define sqr(x) ((x)*(x))
+
/* Implemented in rspl.c: */
extern void alloc_grid(rspl *s);
***************
*** 229,234 ****
--- 237,587 ----
static void solve_gres(mgtmp *m, double tol, int final);
static void init_soln(mgtmp *m1, mgtmp *m2);
static void comp_extrafit_err(mgtmp *m);
+ static int fit_rspl_imp(rspl *s, int flags, void *d, int dtp, int dno,
+ ratai glow, ratai ghigh, int gres[MXDI], ratao vlow, ratao vhigh,
+ double smooth, double avgdev, double weak, void *dfctx,
+ void (*dfunc)(void *cbntx, double *out, double *in));
+
+ /*
+ * Crude optimization of the smoothing factor by minimizing GCV
+ */
+
+ static int
+ opt_gcv(
+ rspl *s, /* this */
+ int flags, /* Combination of flags */
+ void *d, /* Array holding position and function values of data
points */
+ int dtp, /* Flag indicating data type, 0 = (co *), 1 = (cow *) */
+ 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 */
+ 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 */
+ double smooth, /* Smoothing factor, nominal = 1.0 */
+ /* (if -ve, overides optimised smoothing, and sets raw
smoothing */
+ /* typically between 1e-7 .. 1e-1) */
+ double avgdev, /* Average Deviation of input values as proportion of
input range, */
+ /* typical value 0.005 ? (aprox. = 0.564 times the
standard deviation) */
+ double weak, /* Weak weighting, nominal = 1.0 */
+ void *dfctx, /* Opaque weak default function context */
+ void (*dfunc)(void *cbntx, double *out, double *in) /*
Function to set from, NULL if none. */
+ ) {
+ int di = s->di;
+ int fdi = s->fdi;
+ int f, sm, up = 0;
+ double best_smth;
+ double best_gcv_smth;
+ struct rspl_gcv best_gcv;
+ double best_gcv_perchan[MXDO];
+ double best_gcv_smth_perchan[MXDO];
+
+ flags &= ~RSPL_OPT_GCV; /* avoid endless recursion */
+ flags |= RSPL_GCV; /* compute GCV */
+
+ /* XXX fixme: is this assumption correct? */
+ if (flags & RSPL_INCREMENTAL) {
+ error("rspl: Cannot optimize smoothness in incremental mode");
+ }
+
+ if (flags & RSPL_VERBOSE) {
+ printf("Estimating smoothing factor which minimizes GCV\n");
+ }
+
+ best_gcv.gcv_sum = 1e9;
+
+ for (f = 0; f < fdi; f++) {
+ best_gcv_smth_perchan[f] = 0.0;
+ best_gcv_perchan[f] = 1e9;
+ }
+
+ /* Heuristic: terminate if criterion starts to increase again. */
+ /* Wait for 2 consecutive up steps, to avoid local minima. */
+
+ for (sm = 6; sm >= -6 && up < 2; sm--)
+ {
+ rspl *s2;
+ double smth = pow(10, sm/2.0); /* 2 steps per decade */
+
+ if (flags & RSPL_VERBOSE) {
+ printf("Trying %g ...\n", smth);
+ }
+
+ if ((s2 = new_rspl(di, fdi)) == NULL) {
+ error("rspl: opt_gcv: Creation of temporary rspl
failed");
+ }
+
+ /* fit scattered data */
+ fit_rspl_imp(s2, flags, d, dtp, dno, glow, ghigh, gres,
+ vlow, vhigh, smth, avgdev, weak, dfctx, dfunc);
+
+ up++;
+
+ /* did we find a new optimum ? */
+
+ if (s2->gcv.gcv_sum < best_gcv.gcv_sum) {
+ best_gcv_smth = smth;
+ best_gcv = s2->gcv;
+ up = 0;
+ }
+
+ for (f = 0; f < fdi; f++) {
+ if (s2->gcv.gcv[f] < best_gcv_perchan[f]) {
+ best_gcv_smth_perchan[f] = smth;
+ best_gcv_perchan[f] = s2->gcv.gcv[f];
+ up = 0;
+ }
+ }
+
+ s2->del(s2);
+ }
+
+ best_smth = best_gcv_smth;
+
+ if (flags & RSPL_VERBOSE)
+ {
+ printf("Smoothing factor which minimizes GCV: %g ",
best_gcv_smth);
+ for (f = 0; f < fdi; f++)
+ printf("%c%g", f==0 ? '[' : ' ',
best_gcv_smth_perchan[f]);
+ printf("]\n");
+ printf(" GCV: var=%g, std=%g\n", best_gcv.gcv_sum,
sqrt(best_gcv.gcv_sum));
+
+ printf("Using GCV optimized smoothness (%g)\n", best_smth);
+
+ if (smooth != 1.0) {
+ printf("Adjusting by the supplied smoothness
correction: %g\n", smooth);
+ printf("Using smoothing factor: %g\n", best_smth *
smooth);
+ }
+ }
+
+ /*
+ * now fit scattered data with optimal smoothness,
+ * adjusted (multiplied) by the given smoothness.
+ */
+
+ return fit_rspl_imp(s, flags, d, dtp, dno, glow, ghigh, gres, vlow,
+ vhigh, best_smth * smooth, avgdev, weak, dfctx,
dfunc);
+ }
+
+ /*
+ * Compute the randomized GCV estimate, based on a
+ * stochastic estimator of the trace of the influence
+ * operator of the regression model.
+ *
+ * This is based on
+ * Wahba, G., Johnson, D. R., Gao, F. and Gong, J. "Adaptive tuning
+ * of numerical weather prediction models: Part I: randomized GCV and
+ * related methods in three and four dimensional data assimilation."
+ * TR 920, April 1994.
+ */
+
+ static void
+ compute_gcv(
+ rspl *s, /* this */
+ int flags, /* Combination of flags */
+ void *d, /* Array holding position and function values of data
points */
+ int dtp, /* Flag indicating data type, 0 = (co *), 1 = (cow *) */
+ 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 */
+ 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 */
+ double smooth, /* Smoothing factor, nominal = 1.0 */
+ /* (if -ve, overides optimised smoothing, and sets raw
smoothing */
+ /* typically between 1e-7 .. 1e-1) */
+ double avgdev, /* Average Deviation of input values as proportion of
input range, */
+ /* typical value 0.005 ? (aprox. = 0.564 times the
standard deviation) */
+ double weak, /* Weak weighting, nominal = 1.0 */
+ void *dfctx, /* Opaque weak default function context */
+ void (*dfunc)(void *cbntx, double *out, double *in) /* Function to
set from, NULL if none. */
+ ) {
+
+ int i, e, f, iter, niter;
+ int di = s->di;
+ int fdi = s->fdi;
+ rspl *s2 = NULL; /* model fitted to perturbed data */
+ void *d2 = NULL; /* perturbed data */
+ double rss_sum;
+ double tr[MXDO]; /* trace of influence matrix */
+ double rv_std[MXDO]; /* stddev of rv for each ouput channel
*/
+
+ static double *rv[MXDO]; /* random vector for perturbing data */
+ static int nrv[MXDO]; /* number of elements in each rv[i] */
+
+ /* XXX fixme: is this assumption correct? */
+ if (s->inc) {
+ for (f = 0; f < fdi; f++) {
+ s->gcv.rss[f] = s->gcv.gcv[f] =
+ s->gcv.effpar[f] = s->gcv.df[f] = 0.0;
+ }
+ s->gcv.gcv_sum = 0.0;
+ printf("\nrspl: Cannot compute GCV in incremental mode\n");
+ return;
+ }
+
+ /* terminate "***********" */
+ if (s->verbose) printf("\n");
+
+ /* compute no. of iterations */
+
+ for (niter = 1; dno * niter < GCV_MIN_POINTS; niter++)
+ ;
+
+ /*
+ * Create random vector(s) for perturbing the data,
+ * from the binary distribution [-1,1].
+ *
+ * Once allocated, we reuse the _same_ vector(s) for
+ * subsequent GCV computations too, in order to get
+ * more consistent result, if GCV is e.g. used to
+ * find the optimal smoothness factor. Grow the
+ * vector(s) if necessary.
+ */
+
+ for (f = 0; f < fdi; f++) {
+ if (rv[f] == NULL) {
+ /* initial allocation of random vector */
+ if ((rv[f] = malloc(sizeof(double) * dno * niter)) ==
NULL)
+ error("rspl: malloc failed -
compute_gcv:rv[][]");
+ for (i = 0; i < dno * niter; i++)
+ rv[f][i] = rand32(0) & 0x80000000 ? 1.0 : -1.0;
+ nrv[f] = dno * niter;
+ }
+ else if (nrv[f] < dno * niter) {
+ /* reallocate with a larger size */
+ if ((rv[f] = realloc(rv[f], sizeof(double) * dno *
niter)) == NULL)
+ error("rspl: realloc failed -
compute_gcv:rv[][]");
+ /* fill the newly allocated slots */
+ for (i = nrv[f]; i < dno * niter; i++)
+ rv[f][i] = rand32(0) & 0x80000000 ? 1.0 : -1.0;
+ nrv[f] = dno * niter;
+ }
+ }
+
+ /* compute the amount of perturbation */
+ /* heuristically I've choosen 0.1% of full scale */
+
+ for (f = 0; f < fdi; f++) {
+ rv_std[f] = 0.001 * fabs(s->d.vw[f]);
+ #ifdef DEBUG
+ printf("rv_std[%d] = %g\n", f, rv_std[f]);
+ #endif
+ }
+
+ /* initialize tr[] and rss[] */
+
+ for (f = 0; f < fdi; f++) {
+ tr[f] = s->gcv.rss[f] = 0.0;
+ }
+
+ for (iter = 0; iter < niter; iter++)
+ {
+ /* perturb data */
+
+ if (dtp == 0) {
+ co *dd2;
+ if ((d2 = dd2 = malloc(sizeof(co) * dno)) == NULL)
+ error("rspl: malloc failed -
compute_gcv:dd2[]");
+ for (i = 0; i < dno; i++) {
+ dd2[i] = ((co *) d)[i];
+ for (f = 0; f < fdi; f++)
+ dd2[i].v[f] += rv_std[f] *
rv[f][i+dno*iter];
+ }
+ }
+ else {
+ cow *dd2;
+ if ((d2 = dd2 = malloc(sizeof(cow) * dno)) == NULL)
+ error("rspl: malloc failed -
compute_gcv:dd2[]");
+ for (i = 0; i < dno; i++) {
+ dd2[i] = ((cow *) d)[i];
+ for (f = 0; f < fdi; f++)
+ dd2[i].v[f] += rv_std[f] *
rv[f][i+dno*iter];
+ }
+ }
+
+ /* fit perturbed data */
+
+ s2 = new_rspl(di, fdi);
+ fit_rspl_imp(s2, flags&(~RSPL_GCV), d2, dtp, dno, glow, ghigh,
+ gres, vlow, vhigh, smooth, avgdev, weak, dfctx,
dfunc);
+
+ /* terminate "***********" */
+ if (s->verbose) printf("\n");
+
+ /* compute RSS and randomized trace estimate */
+
+ for (i = 0; i < dno; i++) {
+ co point, point1, point2;
+
+ if (dtp == 0) {
+ point1 = ((co *) d)[i];
+ point2 = ((co *) d2)[i];
+ }
+ else {
+ for (e = 0; e < di; e++) {
+ point1.p[e] = ((cow *) d)[i].p[e];
+ point2.p[e] = ((cow *) d2)[i].p[e];
+ }
+ for (f = 0; f < fdi; f++) {
+ point1.v[f] = ((cow *) d)[i].v[f];
+ point2.v[f] = ((cow *) d2)[i].v[f];
+ }
+ }
+
+ point = point1;
+
+ s->interp(s, &point1);
+ s2->interp(s2, &point2);
+
+ /* compute randomized trace estimate */
+
+ for (f = 0; f < fdi; f++) {
+ tr[f] += (1.0/rv_std[f]) * rv[f][i+dno*iter] *
+ (rv_std[f] * rv[f][i+dno*iter] -
(point2.v[f] - point1.v[f]));
+ if (iter == 0) {
+ s->gcv.rss[f] += sqr(point1.v[f] -
point.v[f]);
+ }
+ }
+ }
+
+ s2->del(s2);
+ }
+
+ /* compute GCV, degrees of freedom, effective # of params */
+
+ s->gcv.gcv_sum = 0.0;
+
+ rss_sum = 0.0;
+ for (f = 0; f < fdi; f++) {
+ rss_sum += s->gcv.rss[f];
+ }
+
+ for (f = 0; f < fdi; f++)
+ {
+ tr[f] /= niter;
+ /* limit to > 0, to avoid division by zero later */
+ if (tr[f] < 0.001) tr[f] = 0.001;
+
+ s->gcv.df[f] = tr[f]; /* degrees of freedom */
+ s->gcv.effpar[f] = dno - tr[f]; /* eff. no. parameters
*/
+
+ s->gcv.gcv[f] = s->gcv.rss[f] * dno / sqr(tr[f]);
+ s->gcv.gcv_sum += s->gcv.gcv[f];
+ }
+
+ if (s->verbose) {
+ for (f = 0; f < fdi; f++)
+ printf("Degrees of freedom [%d]: %g\n", f,
s->gcv.df[f]);
+ for (f = 0; f < fdi; f++)
+ printf("Effective number of paramters [%d]: %g\n", f,
s->gcv.effpar[f]);
+ for (f = 0; f < fdi; f++)
+ printf("GCV[%d]: var=%g, std=%g\n", f, s->gcv.gcv[f],
sqrt(s->gcv.gcv[f]));
+ printf("GCV: var=%g, std=%g\n", s->gcv.gcv_sum,
sqrt(s->gcv.gcv_sum));
+ }
+
+ if (d2) free(d2);
+ }
/* Initialise the regular spline from scattered data */
/* Return non-zero if non-monotonic */
***************
*** 256,261 ****
--- 609,615 ----
int di = s->di, fdi = s->fdi;
int i, n, e, f;
int nigc;
+ int retval;
#if defined(__IBMC__) && defined(_M_IX86)
_control87(EM_UNDERFLOW, EM_UNDERFLOW);
***************
*** 267,272 ****
--- 621,632 ----
if (fdi > MXRO)
error("rspl: fit can't handle fdi = %d",fdi);
+ /* optimize smoothing */
+ if (flags & RSPL_OPT_GCV) {
+ return opt_gcv(s, flags, d, dtp, dno, glow, ghigh, gres, vlow,
+ vhigh, smooth, avgdev, weak, dfctx, dfunc);
+ }
+
/* set debug level */
s->debug = (flags >> 24);
***************
*** 469,475 ****
}
/* Do the data point fitting */
! return add_rspl_imp(s, 0, d, dtp, dno);
}
/* Do the work of initialising from initial or extra data points. */
--- 829,843 ----
}
/* Do the data point fitting */
! retval = add_rspl_imp(s, 0, d, dtp, dno);
!
! /* Compute GCV, if requested */
! if (flags & RSPL_GCV) {
! compute_gcv(s, flags, d, dtp, dno, glow, ghigh, gres, vlow,
! vhigh, smooth, avgdev, weak, dfctx, dfunc);
! }
!
! return retval;
}
/* Do the work of initialising from initial or extra data points. */
*** xicc/xicc.h.orig 2006-05-26 06:06:27.000000000 +0200
--- xicc/xicc.h 2006-06-23 20:56:45.000000000 +0200
***************
*** 203,208 ****
--- 203,210 ----
#define ICX_SET_BLACK 0x0002 /* find and set the
black point */
#define ICX_NO_IN_LUTS 0x0040 /* If LuLut: Don't
create input (Device) curves. */
#define ICX_NO_OUT_LUTS 0x0080 /* If LuLut: Don't
create output (PCS) curves. */
+ #define ICX_GCV 0x0100 /* Compute GCV */
+ #define ICX_OPT_GCV 0x0200 /* Optimze smoothness
for minimal GCV */
#define ICX_EXTRA_FIT 0x0400 /* If LuLut: Don't
create output (PCS) curves. */
#define ICX_VERBOSE 0x8000 /* Turn on verboseness
during creation */
struct _icxLuBase * (*set_luobj) (struct _xicc *p,
*** xicc/xlut.c.orig 2006-05-26 06:06:27.000000000 +0200
--- xicc/xlut.c 2006-06-23 22:57:36.000000000 +0200
***************
*** 29,35 ****
#undef DEBUG /* Verbose debug information */
#undef DEBUG_PLOT /* Plot in & out curves */
! #undef INK_LIMIT_TEST /* Turn off input tables for ink limit testing */
#undef CHECK_ILIMIT /* Do sanity checks on meeting ink limit */
#undef NODDV /* Use slow non d/dv powell */
--- 29,35 ----
#undef DEBUG /* Verbose debug information */
#undef DEBUG_PLOT /* Plot in & out curves */
! #undef INK_LIMIT_TEST /* Turn off input tables for ink limit testing
*/
#undef CHECK_ILIMIT /* Do sanity checks on meeting ink limit */
#undef NODDV /* Use slow non d/dv powell */
***************
*** 42,49 ****
#define SHAPE_BASE 0.00001 /* 0 & 1 harmonic weight */
#define SHAPE_HBASE 0.0002 /* 2nd and higher additional weight */
! #undef HACK_ILORD /* 10 */ /* Override input per channel
curve order */
! #undef HACK_OLORD /* 6 */ /* Override output per channel
curve order */
#undef HACK_ISHAPE /* 0.0 */ /* Override input curve
shape */
--- 42,49 ----
#define SHAPE_BASE 0.00001 /* 0 & 1 harmonic weight */
#define SHAPE_HBASE 0.0002 /* 2nd and higher additional weight */
! #undef HACK_ILORD /* 10 */ /* Override input per channel curve order */
! #undef HACK_OLORD /* 6 */ /* Override output per channel curve order */
#undef HACK_ISHAPE /* 0.0 */ /* Override input curve
shape */
***************
*** 2815,2820 ****
--- 2815,2826 ----
double y1[XRES];
#endif /* DEBUG */
+ if (flags & ICX_GCV)
+ rsplflags |= (RSPL_GCV|RSPL_VERBOSE);
+
+ if (flags & ICX_OPT_GCV)
+ rsplflags |= RSPL_OPT_GCV;
+
if (flags & ICX_VERBOSE)
rsplflags |= RSPL_VERBOSE;
***************
*** 3434,3451 ****
/* Initialise from scattered data */
/* Return non-zero if result is non-monotonic */
/* ~~~~ should warn if non-monotonic ??? ~~~~ */
p->clutTable->fit_rspl(
p->clutTable, /* this */
! rsplflags, /* Combination of flags */
! points, /* Array holding position and
function values of data points */
! nodp, /* Number of data points */
! p->ninmin, /* Grid low scale - will expand
to enclose data, NULL = default 0.0 */
! p->ninmax, /* Grid high scale - will
expand to enclose data, NULL = default 1.0 */
! gres, /* Spline grid resolution,
ncells = gres-1 */
! p->noutmin, /* Data value low normalize,
NULL = default 0.0 */
! p->noutmax, /* Data value high normalize -
NULL = default 1.0 */
! smooth, /* Smoothing factor, nominal =
1.0 */
! avgdev /* reading Average Deviation as
a proportion of the input range */
);
}
--- 3440,3458 ----
/* Initialise from scattered data */
/* Return non-zero if result is non-monotonic */
/* ~~~~ should warn if non-monotonic ??? ~~~~ */
+
p->clutTable->fit_rspl(
p->clutTable, /* this */
! rsplflags, /* Combination of flags */
! points, /* Array holding position and function
values of data points */
! nodp, /* Number of data points */
! p->ninmin, /* Grid low scale - will expand to
enclose data, NULL = default 0.0 */
! p->ninmax, /* Grid high scale - will expand to
enclose data, NULL = default 1.0 */
! gres, /* Spline grid resolution, ncells =
gres-1 */
! p->noutmin, /* Data value low normalize, NULL =
default 0.0 */
! p->noutmax, /* Data value high normalize - NULL =
default 1.0 */
! smooth, /* Smoothing factor, nominal = 1.0 */
! avgdev /* reading Average Deviation as a
proportion of the input range */
);
}
*** profile/prof.h.orig 2006-06-23 23:00:48.000000000 +0200
--- profile/prof.h 2006-06-23 20:56:45.000000000 +0200
***************
*** 68,73 ****
--- 68,74 ----
int fwacomp, /* FWA compensation requested */
double smooth, /* RSPL smoothing factor, -ve if raw */
double avgdev, /* reading Average Deviation as a
proportion of the input range */
+ int opt_cv, /* flag, whether to optimize
smoothness, 1=GCV, others TBD. */
char *ipname, /* input icc profile - enables gamut
map, NULL if none */
char *sgname, /* source image gamut - NULL if none */
char *absname, /* abstract profile name - NULL if none
*/
***************
*** 92,97 ****
--- 93,99 ----
cgats *icg, /* input cgats structure */
double smooth, /* RSPL smoothing factor, -ve if raw */
double avgdev, /* reading Average Deviation as a
proportion of the input range */
+ int opt_cv, /* flag, whether to optimize
smoothness, 1=GCV, others TBD. */
profxinf *pi /* Optional Profile creation extra data
*/
);
*** profile/profile.c.orig 2006-05-26 06:06:35.000000000 +0200
--- profile/profile.c 2006-06-23 21:00:29.000000000 +0200
***************
*** 98,103 ****
--- 98,105 ----
fprintf(stderr," 1931_2, 1964_10, S&B 1955_2, 1964_10c,
shaw, J&V 1978_2 (def.)\n");
fprintf(stderr," -f Use Fluorescent Whitening Agent
compensation\n");
fprintf(stderr," -r avgdev Average deviation of device+instrument
readings as a percentage (default %4.2f%%)\n",DEFAVGDEV);
+ fprintf(stderr," -O opt Optimize smoothing factor:\n");
+ fprintf(stderr," 1 = minimize GCV\n");
/* fprintf(stderr," -rs smooth RSPL suplimental optimised smoothing
factor\n"); */
/* fprintf(stderr," -rr smooth RSPL raw underlying smoothing
factor\n"); */
fprintf(stderr," -s src.icc Apply gamut mapping to perceptual B2A
table for given source space\n");
***************
*** 148,153 ****
--- 150,156 ----
int fwacomp = 0; /* FWA compensation */
double avgdev = DEFAVGDEV/100.0; /* Average measurement
deviation */
double smooth = 1.0; /* RSPL Smoothness factor (relative,
for verification) */
+ int opt_cv = 0; /* optimize smoothness for min.
GCV, others TBD */
int spec = 0; /* Use spectral data flag */
icxIllumeType illum = icxIT_D50; /* Spectral defaults */
xspect cust_illum; /* Custom illumination spectrum
*/
***************
*** 175,180 ****
--- 178,185 ----
error_program = argv[0];
memset((void *)&xpi, 0, sizeof(profxinf)); /* Init extra profile
info to defaults */
+ rand32(time(NULL));
+
/* Init VC overrides so that we know when the've been set */
ivc_p.Ev = -1;
ivc_p.Wxyz[0] = -1.0; ivc_p.Wxyz[1] = -1.0; ivc_p.Wxyz[2] = -1.0;
***************
*** 221,227 ****
usage("Usage requested");
else if (argv[fa][1] == 'v' || argv[fa][1] == 'V')
! verb = 1;
/* Manufacturer description string */
else if (argv[fa][1] == 'A') {
--- 226,232 ----
usage("Usage requested");
else if (argv[fa][1] == 'v' || argv[fa][1] == 'V')
! verb++;
/* Manufacturer description string */
else if (argv[fa][1] == 'A') {
***************
*** 458,464 ****
}
/* Spectral Observer type */
! else if (argv[fa][1] == 'o' || argv[fa][1] == 'O') {
fa = nfa;
if (na == NULL) usage("Expect argument to
observer flag -o");
if (strcmp(na, "1931_2") == 0) {
/* Classic 2 degree */
--- 463,469 ----
}
/* Spectral Observer type */
! else if (argv[fa][1] == 'o') {
fa = nfa;
if (na == NULL) usage("Expect argument to
observer flag -o");
if (strcmp(na, "1931_2") == 0) {
/* Classic 2 degree */
***************
*** 503,508 ****
--- 508,523 ----
}
}
+ /* optimize smoothness */
+ else if (argv[fa][1] == 'O') {
+ fa = nfa;
+ if (na == NULL) usage("Expect argument to
optimization flag -O");
+ opt_cv = atoi(na);
+ if (opt_cv < 1 || opt_cv > 1) {
+ usage("Illegal argument to optimization
flag, -O1 expected");
+ }
+ }
+
/* Percetual Source Gamut and Perceptual/Saturation
Gamut Maping mode enable */
else if (argv[fa][1] == 's' || argv[fa][1] == 'S') {
if (argv[fa][1] == 'S')
***************
*** 785,791 ****
}
make_output_icc(ptype, verb, iquality, oquality, noiluts,
nooluts, verify, &ink, outname,
! icg, spec, illum, &cust_illum, observ, fwacomp,
smooth, avgdev,
ipname[0] != '\000' ? ipname : NULL,
sgname[0] != '\000' ? sgname : NULL,
absname[0] != '\000' ? absname : NULL,
--- 800,806 ----
}
make_output_icc(ptype, verb, iquality, oquality, noiluts,
nooluts, verify, &ink, outname,
! icg, spec, illum, &cust_illum, observ, fwacomp,
smooth, avgdev, opt_cv,
ipname[0] != '\000' ? ipname : NULL,
sgname[0] != '\000' ? sgname : NULL,
absname[0] != '\000' ? absname : NULL,
***************
*** 796,802 ****
if (ptype == prof_default)
ptype = prof_clutLab; /* For best possible
quality */
! make_input_icc(ptype, verb, iquality, verify, nsabs, outname,
icg, smooth, avgdev, &xpi);
} else if (strcmp(icg->t[0].kdata[ti],"DISPLAY") == 0) {
--- 811,817 ----
if (ptype == prof_default)
ptype = prof_clutLab; /* For best possible
quality */
! make_input_icc(ptype, verb, iquality, verify, nsabs, outname,
icg, smooth, avgdev, opt_cv, &xpi);
} else if (strcmp(icg->t[0].kdata[ti],"DISPLAY") == 0) {
***************
*** 809,815 ****
/* ICC V2.3 doesn't have display intents. */
/* ICC V2.4 does. What should we do here ? */
make_output_icc(ptype, verb, iquality, oquality, noiluts,
nooluts, verify, NULL, outname,
! icg, spec, illum, &cust_illum, observ, 0,
smooth, avgdev,
NULL, NULL,
absname[0] != '\000' ? absname : NULL,
0, &ivc_p, &ovc_p, ivc_e, ovc_e,
--- 824,830 ----
/* ICC V2.3 doesn't have display intents. */
/* ICC V2.4 does. What should we do here ? */
make_output_icc(ptype, verb, iquality, oquality, noiluts,
nooluts, verify, NULL, outname,
! icg, spec, illum, &cust_illum, observ, 0,
smooth, avgdev, opt_cv,
NULL, NULL,
absname[0] != '\000' ? absname : NULL,
0, &ivc_p, &ovc_p, ivc_e, ovc_e,
*** profile/profin.c.orig 2006-06-23 23:01:00.000000000 +0200
--- profile/profin.c 2006-06-23 20:56:45.000000000 +0200
***************
*** 81,86 ****
--- 81,87 ----
cgats *icg, /* input cgats structure */
double smooth, /* RSPL smoothing factor, -ve if raw */
double avgdev, /* reading Average Deviation as a
proportion of the input range */
+ int opt_cv, /* flag, whether to optimize
smoothness, 1=GCV, others TBD. */
profxinf *xpi /* Optional Profile creation extra data
*/
) {
icmFile *wr_fp;
***************
*** 445,453 ****
--- 446,460 ----
if ((wr_xicc = new_xicc(wr_icco)) == NULL)
error ("Creation of xicc failed");
+ if (opt_cv == 1)
+ flags |= ICX_OPT_GCV;
+
if (verb)
flags |= ICX_VERBOSE;
+ if (verb >= 2)
+ flags |= ICX_GCV;
+
if (nsabs == 0)
flags |= ICX_SET_WHITE | ICX_SET_BLACK; /* Set white
and black */
*** profile/profout.c.orig 2006-05-26 06:06:35.000000000 +0200
--- profile/profout.c 2006-06-23 20:56:45.000000000 +0200
***************
*** 446,451 ****
--- 446,452 ----
int fwacomp, /* FWA compensation requested */
double smooth, /* RSPL smoothing factor, -ve if raw */
double avgdev, /* reading Average Deviation as a
proportion of the input range */
+ int opt_cv, /* flag, whether to optimize
smoothness, 1=GCV, others TBD. */
char *ipname, /* input icc profile - enables gamut
map, NULL if none */
char *sgname, /* source image gamut - NULL if none */
char *absname, /* abstract profile name - NULL if none
*/
***************
*** 1415,1423 ****
--- 1416,1430 ----
if (nooluts)
flags |= ICX_NO_OUT_LUTS;
+ if (opt_cv == 1)
+ flags |= ICX_OPT_GCV;
+
if (verb)
flags |= ICX_VERBOSE;
+ if (verb >= 2)
+ flags |= ICX_GCV;
+
/* Setup Device -> PCS conversion (Fwd) object from
scattered data. */
if ((AtoB = wr_xicc->set_luobj(
wr_xicc, icmFwd, isdisp ?
icmDefaultIntent : icRelativeColorimetric,
*** spectro/fakeread.c.orig 2006-05-26 06:06:31.000000000 +0200
--- spectro/fakeread.c 2006-06-23 20:56:45.000000000 +0200
***************
*** 813,820 ****
double dv = PCS[j];
double rr = d_rand(-0.5 * rplevel, 0.5
* rplevel);
dv += rr;
! if (dv < 0.0)
dv = 0.0;
PCS[j] = dv;
}
}
--- 813,823 ----
double dv = PCS[j];
double rr = d_rand(-0.5 * rplevel, 0.5
* rplevel);
dv += rr;
!
! /* don't enforce b* or a* >= 0 !!! */
! if (dv < 0.0 && (!dolab || j == 0))
dv = 0.0;
+
PCS[j] = dv;
}
}
Other related posts:
