gronod/argyllcms
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 11 Wiki Activity
Files
development
argyllcms/xicc/xfit.c
T
gronod 45bdf51201 Initial commit of v3.5.0
2026-08-20 20:28:40 +01:00

3151 lines
90 KiB
C
Raw Permalink Blame History

/*
* International Color Consortium color transform expanded support
*
* Author: Graeme W. Gill
* Date: 2/7/00
* Version: 1.00
*
* Copyright 2000 - 2007 Graeme W. Gill
* All rights reserved.
* This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
* see the License.txt file for licencing details.
*
* Based on the xlut.c code.
*/
/*
* TTBD:
*
#ifdef PRIVATE
* [ This is probably key to fixing L*a*b* B2A black bumpiness ??. ]
*
* Need to create something similar to this for B2A tables, device link
* creation and imdi creation, rather than using the existing per
* channel curves or their inverses. Call it pchfit ?
* Would also be nice to add grid sizing
* to cover just gamut range ? - just spec. it as parameter.
*
* !! This might eliminate any need for dividing transforms into
* per-channel/core/per-channel parts in icc !!
*
* Might work just as well as this xfit code if we create prelim cLUT
* and use this new routine to create input & output values ? (worth testing, but not essential).
*
* Idea would be to treat target transform as a black box, plus any needed
* conversions from device to PCS.
*
* For the input the aim would be to distribute cLut index points
* perceptually evenly - i.e. like the print calibrate code.
* For Lab in this is a NOP, for XYZ in this would be an XYZ to L* curve,
* and for device in this would be a dev->delta E curve. Possibly:
* Run set of device input values through ->Lab, and then create
* smart non-mon curve fit. Make it smart by first fitting
* power curve from black end vs. normal curve, picking best
* and then fitting higher orders from that choice (or fix curves to be auto smart
* about fitting curves ?)
*
* For the output curve the aim would be to create a non-extreme
* mapping that makes the output space more linear mixing.
* One approach would be to start with a default curve that maps
* output to a linear light space, i.e. XYZ would be a NOP,
* Lab would be a L* to Y curve, device would be device to XYZ curves,
* then add mod curve to improve additive behavior.
*
* Optimization target would be sampling points along the neutral +
* primary axes. Use adaptive point sampling to arrive at (say) 50-100 point
* along each axis roughly evenly spaced.
* Or take cLUT grid res as parameter, and make test points center points
* of critical color vectors cells ?
* The smart non-mon curves
* that minimizes interleave linear interp errors between points on each axis.
* Need to make a*b* space symetrical about zero.
* Also need to treat non diagonal device spaces specially ?
*
#endif
*
* Currently the Lab A2B output tables are adjusted for ab symetry
* to make the B2A white point land on a grid point, given that
* the icc code forces a symetric ab range. This only works
* because the B2A is using the same per channel curves.
* Done properly, it should be possible to know where grid
* points land within the range, and to modify the B2A input curves
* to make the white point land on a grid point.
* (For B2A the pseudo least squares adjustment needs to be turned
* off for that grid point too.)
*
* Note that one quandry is that the curve fitting doesn't
* fit well when the input data has an offset and/or plateaus.
*/
/*
* This module provides curve and matrix fitting functionality used to
* create per channel input and/or output curves for clut profiles,
* and optionally creates the rspl to fit the data as well.
*
* The approach used is to initialy creating input and output shaper
* curves that minimize the overall delta E of the test point set,
* when a linear matrix is substituted for a clut is used.
* (This is the same as pre-V0.70 approach. )
*
* The residual error between the test point set and the shaper/matrix/shaper
* model is then computed, and mapped against each input channel, and
* a positioning mapping curve created that aims to map rspl grid
* locations more densly where residual errors are high, and more
* sparsely where they are low.
* The input shaper and positioning curves are then combined together,
* so that the positioning curve determines which input values map
* to rspl grid points, and the shaper curve determines the mapping
* between the grid points. The width of the grid cells is computed
* in the shaper mapped input space, and then fed to the rspl fitting
* code, so that its smoothness evaluation function can take account
* of the non-uniform spacing of the grid points.
*/
#include <sys/types.h>
#include <string.h>
#include <ctype.h>
#ifdef __sun
#include <unistd.h>
#endif
#include "copyright.h"
#include "aconfig.h"
#include "numlib.h"
#include "icc.h"
#include "rspl.h"
#include "plot.h"
#include "conv.h"
#include "xicc.h"
#include "xfit.h"
#include "sort.h"
#undef USE_XYZ_Y2LCURVE /* [Und] Use underlying L* curve for XYZ encoding */
/* This seems to work badly, even with high smoothness. Why ? */
/* It does speed up 1D lut creation though. */
#undef DEBUG /* Debug information */
#undef DEBUG_PROGRESS /* Show powell progress */
#undef DEBUG_PLOT /* Plot in & out curves */
#undef SPECIAL_FORCE /* Check rspl nodes against linear XYZ model */
#undef SPECIAL_FORCE_GAMMA /* Force correct gamma shaper curves */
#undef SPECIAL_TEST_GAMMA /* Use gamma on reference model */
#undef SPECIAL_TEST_LAB
#define RSPLFLAGS (0 /* | RSPL_2PASSSMTH */ /* | RSPL_EXTRAFIT2 */)
#undef EXTEND_GRID /* [Undef] Use extended RSPL grid around interpolation */
#define EXTEND_GRID_BYN 2 /* Rows to extend grid. */
#undef NODDV /* Use slow non d/dv powell else use conjgrad */
#define CURVEPOW 1.0 /* Power to raise deltaE squared to in setting in/out curves */
/* This provides a means of punishing high maximum errors. */
#define POWTOL1 1e-3 /* Shaper Powell optimiser tollerance for first passes */
#define MAXITS1 1000 /* Shaper number of itterations for first passes */
#define POWTOL 1e-5 /* Shaper Powell optimiser tollerance in delta E squared ^ CURVEPOW */
#define MAXITS 4000 /* Shaper number of itterations before giving up */
#define PDDEL 1e-6 /* Fake partial derivative del */
/* Weights for shaper in/out curve parameters, to minimise unconstrained "wiggles" */
#define SHAPE_WEIGHT 1.0 /* Overal shaper weight contribution - err on side of smoothness */
#define SHAPE_HW01 0.002 /* 0 & 1 harmonic weights */
#define SHAPE_HBREAK 4 /* Harmonic that has HWBR */
#define SHAPE_HWBR 20.0 /* Base weight of harmonics HBREAK up */
#define SHAPE_HWINC 60.0 /* Increase in weight for each harmonic above HWBR */
/* Weights for the positioning curve parameters */
#define PSHAPE_MINE 0.02 /* Minum background residual error level */
#define PSHAPE_DIST 1.0 /* Agressivness of grid distribution */
/* - - - - - - - - - - - - - - - - - */
/* Extra non-linearity used as base for XYZ output curves. */
/* This makes the XYZ grid values more perceptual, and asks less */
/* of the automatically created output curve shape. */
/* (We assume XYZ is in 0..1 scale */
/* Transfer function with offset and scale + Y2L curve */
static double icxSTransFuncY2L(
double *v, /* Pointer to first parameter */
int luord, /* Number of parameters */
double vv, /* Source of value */
double min, /* Scale values */
double max
) {
max -= min;
vv = (vv - min)/max;
#ifdef USE_XYZ_Y2LCURVE
if (vv > 0.008856451586)
vv = 1.16 * pow(vv,1.0/3.0) - 0.16;
else
vv = 9.032962896 * vv;
#endif
vv = icxTransFunc(v, luord, vv);
vv = (vv * max) + min;
return vv;
}
/* Inverse Transfer function with offset and scale + Y2L */
static double icxInvSTransFuncY2L(
double *v, /* Pointer to first parameter */
int luord, /* Number of parameters */
double vv, /* Source of value */
double min, /* Scale values */
double max
) {
max -= min;
vv = (vv - min)/max;
vv = icxInvTransFunc(v, luord, vv);
#ifdef USE_XYZ_Y2LCURVE
if (vv > 0.08)
vv = pow((vv + 0.16)/1.16, 3.0);
else
vv = vv/9.032962896;
#endif
vv = (vv * max) + min;
return vv;
}
/* Transfer function with offset and scale, and */
/* partial derivative with respect to the */
/* parameters and the input value. */
static double icxdpdiSTransFuncY2L(
double *v, /* Pointer to first parameter */
double *dv, /* Return derivative wrt each parameter */
double *pdin, /* Return derivative wrt source value */
int luord, /* Number of parameters */
double vv, /* Source of value */
double min, /* Scale values */
double max
) {
int i;
double idv = 1.0;
max -= min;
#ifdef USE_XYZ_Y2LCURVE
if (vv > 0.008856451586) {
vv = 1.16 * pow(vv,1.0/3.0) - 0.16;
idv = 1.16 / (3.0 * pow(vv, 2.0/3.0));
} else {
vv = 9.032962896 * vv;
idv = 9.032962896;
}
#endif
vv = (vv - min)/max;
vv = icxdpdiTransFunc(v, dv, pdin, luord, vv);
*pdin *= idv; /* Account for input multiplier */
vv = (vv * max) + min;
for (i = 0; i < luord; i++) {
dv[i] *= max;
}
return vv;
}
/* - - - - - - - - - - - - - - - - - */
#ifdef DEBUG
static void dump_xfit(xfit *p);
#endif
#ifdef DEBUG_PLOT /* Not currently used in runtime code*/
/* Lookup a value though an input position curve */
static double xfit_poscurve(xfit *p, double in, int chan) {
double rv = in;
if (p->tcomb & oc_p)
rv = icxSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], rv,
p->in_min[chan], p->in_max[chan]);
return rv;
}
/* Lookup di values though the input position curves */
static void xfit_poscurves(xfit *p, double *out, double *in) {
int e;
for (e = 0; e < p->di; e++) {
double val = in[e];
if (p->tcomb & oc_i)
val = icxSTransFunc(p->v + p->pos_offs[e], p->iluord[e], val,
p->in_min[e], p->in_max[e]);
out[e] = val;
}
}
/* Inverse Lookup a value though an input position curve */
static double xfit_invposcurve(xfit *p, double in, int chan) {
double rv = in;
if (p->tcomb & oc_i)
rv = icxInvSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], rv,
p->in_min[chan], p->in_max[chan]);
return rv;
}
/* Inverse Lookup di values though input position curves */
static void xfit_invposcurves(xfit *p, double *out, double *in) {
int e;
for (e = 0; e < p->di; e++) {
double val = in[e];
if (p->tcomb & oc_i)
val = icxInvSTransFunc(p->v + p->pos_offs[e], p->iluord[e], val,
p->in_min[e], p->in_max[e]);
out[e] = val;
}
}
#endif /* DEBUG_PLOT */
/* - - - - - - - - - - - - - - - - - */
/* Lookup a value though input shape curve */
static double xfit_shpcurve(xfit *p, double in, int chan) {
double rv = in;
if (p->tcomb & oc_i) {
#ifdef SPECIAL_FORCE_GAMMA
double gam;
if (chan == 0)
gam = 1.9;
else if (chan == 1)
gam = 2.0;
else if (chan == 2)
gam = 2.1;
else
gam = 1.0;
rv = pow(in, gam);
#else
rv = icxSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], rv,
p->in_min[chan], p->in_max[chan]);
#endif
}
return rv;
}
#ifdef NEVER /* Not currently used */
/* Lookup a value though shape curves */
static void xfit_shpcurves(xfit *p, double *out, double *in) {
int e;
for (e = 0; e < p->di; e++) {
double val = in[e];
if (p->tcomb & oc_i)
val = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], val,
p->in_min[e], p->in_max[e]);
out[e] = val;
}
}
#endif /* NEVER */
/* Inverse Lookup a value though a shape curve */
static double xfit_invshpcurve(xfit *p, double in, int chan) {
double rv = in;
if (p->tcomb & oc_i) {
#ifdef SPECIAL_FORCE_GAMMA
double gam;
if (chan == 0)
gam = 1.9;
else if (chan == 1)
gam = 2.0;
else if (chan == 2)
gam = 2.1;
else
gam = 1.0;
rv = pow(rv, 1.0/gam);
#else
rv = icxInvSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], rv,
p->in_min[chan], p->in_max[chan]);
#endif
}
return rv;
}
#ifdef NEVER /* Not currently used */
/* Inverse Lookup a value though shape curves */
static void xfit_invshpcurves(xfit *p, double *out, double *in) {
int e;
for (e = 0; e < p->di; e++) {
double val = in[e];
if (p->tcomb & oc_i)
val = icxInvSTransFunc(p->v + p->shp_offs[e], p->iluord[e], val,
p->in_min[e], p->in_max[e]);
out[e] = val;
}
}
#endif /* NEVER */
/* - - - - - - - - - - - - - - - - - */
/* Lookup values through the shaper/matrix/shaper model */
static void xfit_shmatsh(xfit *p, double *out, double *in) {
double tin[MXDI];
int e, f;
for (e = 0; e < p->di; e++)
tin[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], in[e],
p->in_min[e], p->in_max[e]);
icxCubeInterp(p->v + p->mat_off, p->fdi, p->di, out, tin);
if (p->flags & XFIT_OUT_LAB) {
for (f = 0; f < p->fdi; f++)
out[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], out[f],
p->out_min[f], p->out_max[f]);
} else {
for (f = 0; f < p->fdi; f++)
out[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], out[f],
p->out_min[f], p->out_max[f]);
}
}
/* - - - - - - - - - - - - - - - - - - - - */
/* Combined input positioning & shaper transfer curve functions */
//int db = 0;
/* Lookup a value though the input positioning and shaper curves */
static double xfit_inpscurve(xfit *p, double in, int chan) {
double rv;
/* Just shaper curve */
if ((p->tcomb & oc_ip) == oc_i) {
rv = icxSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], in,
p->in_min[chan], p->in_max[chan]);
/* shaper and positioning */
} else if ((p->tcomb & oc_ip) == oc_ip) {
double nin, npind; /* normalized in, normalized position in' */
int six; /* Span index */
double npind0, npind1; /* normalized position in' span values */
double npin0, npin1; /* normalized position in span values */
double nsind0, nsind1; /* normalized shaper in' span values */
double nsind; /* normalised shaper in' value */
/* Normalize */
nin = (in - p->in_min[chan])/(p->in_max[chan] - p->in_min[chan]);
//if (db) printf("\n~1 inpscurve: cha %d, input value %f, norm %f\n",chan,in,nin);
/* Locate the span the input point will be in after positioning lookup */
npind = icxTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], nin);
/* Quantize position space value to grid */
six = (int)floor(npind * (p->gres[chan]-1.0));
if (six > (p->gres[chan]-2))
six = (p->gres[chan]-2);
/* Compute span position position in' values */
npind0 = six / (p->gres[chan]-1.0);
npind1 = (six + 1.0) / (p->gres[chan]-1.0);
//if (db) printf("~1 npind %f, six %d, npind0 %f, npind1 %f\n",npind,six,npind0,npind1);
/* Compute span in values */
npin0 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind0);
npin1 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind1);
//if (db) printf("~1 npin0 %f, npin1 %f\n",npin0,npin1);
/* Compute shaper space values of in' and spane */
#ifdef NEVER
nsind = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], nin);
nsind0 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin0);
nsind1 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin1);
#else
nsind = xfit_shpcurve(p, nin, chan);
nsind0 = xfit_shpcurve(p, npin0, chan);
nsind1 = xfit_shpcurve(p, npin1, chan);
#endif
//if (db) printf("~1 nsind %f, nsind0 %f, nsind1 %f\n",nsind,nsind0,nsind1);
/* Offset and scale shaper in' value to match position span */
rv = (nsind - nsind0)/(nsind1 - nsind0) * (npind1 - npind0) + npind0;
//if (db) printf("~1 scale offset ind %f\n",rv);
/* de-normalize */
rv = rv * (p->in_max[chan] - p->in_min[chan]) + p->in_min[chan];
//if (db) printf("~1 returning %d\n",rv);
/* Just positioning curve */
} else if ((p->tcomb & oc_ip) == oc_p) {
rv = icxSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], in,
p->in_min[chan], p->in_max[chan]);
} else {
rv = in;
}
return rv;
}
/* Lookup di values though the input positioning and shaper curves */
static void xfit_inpscurves(xfit *p, double *out, double *in) {
int e;
for (e = 0; e < p->di; e++)
out[e] = xfit_inpscurve(p, in[e], e);
}
/* Inverse Lookup a value though the input positioning and shaper curves */
static double xfit_invinpscurve(xfit *p, double in, int chan) {
double rv;
/* Just shaper curve */
if ((p->tcomb & oc_ip) == oc_i) {
rv = icxInvSTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], in,
p->in_min[chan], p->in_max[chan]);
/* shaper and positioning */
} else if ((p->tcomb & oc_ip) == oc_ip) {
double nind, nin; /* normalized in', normalized in */
int six; /* Span index */
double npind0, npind1; /* normalized position in' span values */
double npin0, npin1; /* normalized position in span values */
double nsind0, nsind1; /* normalized shaper in' span values */
double nsind; /* normalized shaper in' value */
/* Normalize */
nind = (in - p->in_min[chan])/(p->in_max[chan] - p->in_min[chan]);
//if (db) printf("\n~1 invinpscurve: cha %d, input value %f, norm %f\n",chan,in,nind);
/* Quantize to grid */
six = (int)floor(nind * (p->gres[chan]-1.0));
if (six > (p->gres[chan]-2))
six = (p->gres[chan]-2);
/* Compute span in' values */
npind0 = six / (p->gres[chan]-1.0);
npind1 = (six + 1.0) / (p->gres[chan]-1.0);
//if (db) printf("~1 six %d, npind0 %f, npind1 %f\n",six,npind0,npind1);
/* Lookup span in values through position curve */
npin0 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind0);
npin1 = icxInvTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], npind1);
//if (db) printf("~1 npin0 %f, npin1 %f\n",npin0,npin1);
/* Compute span shaper in' values */
#ifdef NEVER
nsind0 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin0);
nsind1 = icxTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], npin1);
#else
nsind0 = xfit_shpcurve(p, npin0, chan);
nsind1 = xfit_shpcurve(p, npin1, chan);
#endif
/* Offset and scale position in' value to match shaper span */
nsind = (nind - npind0)/(npind1 - npind0) * (nsind1 - nsind0) + nsind0;
//if (db) printf("~1 nsind %f, nsind0 %f, nsind1 %f\n",nsind,nsind0,nsind1);
/* Invert through shaper curve */
#ifdef NEVER
nin = icxInvTransFunc(p->v + p->shp_offs[chan], p->iluord[chan], nsind);
#else
nin = xfit_invshpcurve(p, nsind, chan);
#endif
/* de-normalize */
rv = nin * (p->in_max[chan] - p->in_min[chan]) + p->in_min[chan];
//if (db) printf("\n~1 nin = %f, returning %f\n",nin,rv);
/* Just positioning curve */
} else if ((p->tcomb & oc_ip) == oc_p) {
rv = icxInvSTransFunc(p->v + p->pos_offs[chan], p->iluord[chan], in,
p->in_min[chan], p->in_max[chan]);
} else {
rv = in;
}
return rv;
}
#ifdef NEVER
/* Check that inverse is working */
static double _xfit_inpscurve(xfit *p, double in, int chan) {
double inv, rv, iinv;
inv = in;
rv = xfit_inpscurve(p, in, chan);
iinv = xfit_invinpscurve(p, rv, chan);
if (fabs(in - iinv) > 1e-5)
warning("xfit_inpscurve check, got %f, should be %f\n",iinv,in);
return rv;
}
#define xfit_inpscurve _xfit_inpscurve
#endif /* NEVER */
/* Inverse Lookup di values though the input positioning and shaper curves */
static void xfit_invinpscurves(xfit *p, double *out, double *in) {
int e;
for (e = 0; e < p->di; e++)
out[e] = xfit_invinpscurve(p, in[e], e);
}
/* - - - - - - - - - - - - - - - - - - - - */
/* Combined output transfer curve functions */
/* Lookup a value though an output curve */
static double xfit_outcurve(xfit *p, double in, int chan) {
double rv;
if (p->tcomb & oc_o) {
if (p->flags & XFIT_OUT_LAB)
rv = icxSTransFunc(p->v + p->out_offs[chan], p->oluord[chan], in,
p->out_min[chan], p->out_max[chan]);
else
rv = icxSTransFuncY2L(p->v + p->out_offs[chan], p->oluord[chan], in,
p->out_min[chan], p->out_max[chan]);
} else
rv = in;
return rv;
}
/* Lookup fdi values though the output curves */
static void xfit_outcurves(xfit *p, double *out, double *in) {
int f;
if (p->flags & XFIT_OUT_LAB) {
for (f = 0; f < p->fdi; f++) {
double val = in[f];
if (p->tcomb & oc_o)
val = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], val,
p->out_min[f], p->out_max[f]);
out[f] = val;
}
} else {
for (f = 0; f < p->fdi; f++) {
double val = in[f];
if (p->tcomb & oc_o)
val = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], val,
p->out_min[f], p->out_max[f]);
out[f] = val;
}
}
}
/* Inverse Lookup a value though an output curve */
static double xfit_invoutcurve(xfit *p, double in, int chan) {
double rv;
if (p->tcomb & oc_o) {
if (p->flags & XFIT_OUT_LAB)
rv = icxInvSTransFunc(p->v + p->out_offs[chan], p->oluord[chan], in,
p->out_min[chan], p->out_max[chan]);
else
rv = icxInvSTransFuncY2L(p->v + p->out_offs[chan], p->oluord[chan], in,
p->out_min[chan], p->out_max[chan]);
} else
rv = in;
return rv;
}
/* Inverse Lookup fdi values though output curves */
static void xfit_invoutcurves(xfit *p, double *out, double *in) {
int f;
if (p->flags & XFIT_OUT_LAB) {
for (f = 0; f < p->fdi; f++) {
double val = in[f];
if (p->tcomb & oc_o)
val = icxInvSTransFunc(p->v + p->out_offs[f], p->oluord[f], val,
p->out_min[f], p->out_max[f]);
out[f] = val;
}
} else {
for (f = 0; f < p->fdi; f++) {
double val = in[f];
if (p->tcomb & oc_o)
val = icxInvSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], val,
p->out_min[f], p->out_max[f]);
out[f] = val;
}
}
}
/* - - - - - - - - - */
/* Convert an output value from absolute */
/* to relative using the current white point. */
static void xfit_abs_to_rel(xfit *p, double *out, double *in) {
if (p->flags & XFIT_OUT_WP_REL) {
if (p->flags & XFIT_OUT_LAB) {
icmLab2XYZ(&icmD50, out, in);
icmMulBy3x3(out, p->fromAbs, out);
icmXYZ2Lab(&icmD50, out, out);
} else {
icmMulBy3x3(out, p->fromAbs, in);
}
} else {
out[0] = in[0];
out[1] = in[1];
out[2] = in[2];
}
}
/* Convert an XYZ output value from absolute */
/* to cLut relative using the current white point. */
static void xfit_XYZ_abs_to_rel(xfit *p, double *out, double *in) {
if (p->flags & XFIT_OUT_WP_REL) {
if (p->flags & XFIT_OUT_LAB) {
icmMulBy3x3(out, p->fromAbs, in);
icmXYZ2Lab(&icmD50, out, out);
} else {
icmMulBy3x3(out, p->fromAbs, in);
}
} else {
if (p->flags & XFIT_OUT_LAB) {
icmXYZ2Lab(&icmD50, out, in);
} else {
out[0] = in[0];
out[1] = in[1];
out[2] = in[2];
}
}
}
/* - - - - - - - - - */
/* return a weighting for the magnitude of the in and out */
/* shaping parameters squared. This is to reduce unconstrained "wiggles" */
static double shapmag(
xfit *p /* Base of optimisation structure */
) {
double tt, w;
double *b; /* Base of parameters for this section */
int di = p->di;
int fdi = p->fdi;
int e, f, k;
double iparam = 0.0;
double oparam = 0.0;
double dd;
if (p->opt_msk & oc_i) {
dd = SHAPE_WEIGHT/(double)(di);
b = p->v + p->shp_off;
for (e = 0; e < di; e++) {
for (k = 0; k < p->iluord[e]; k++) {
if (k <= 1) {
w = SHAPE_HW01;
} else if (k <= SHAPE_HBREAK) {
double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
w *= p->shp_smooth[e];
} else {
w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
w *= p->shp_smooth[e];
}
tt = *b++;
tt *= tt; /* Squared */
iparam += w * tt;
}
}
iparam *= dd;
}
if (p->opt_msk & oc_o) {
dd = SHAPE_WEIGHT/(double)(fdi);
b = p->v + p->out_off;
for (f = 0; f < fdi; f++) {
for (k = 0; k < p->oluord[f]; k++) {
if (k <= 1) {
w = SHAPE_HW01;
} else if (k <= SHAPE_HBREAK) {
double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
w *= p->out_smooth[f];
} else {
w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
w *= p->out_smooth[f];
}
tt = *b++;
tt *= tt; /* Squared */
oparam += w * tt;
}
}
oparam *= dd;
}
return iparam + oparam;
}
/* return a weighting for the magnitude of the in and out */
/* shaping parameters. This is to reduce unconstrained "wiggles" */
/* Also sum the partial derivative for the parameters involved */
static double dshapmag(
xfit *p, /* Base of optimisation structure */
double *dav /* Sum del's */
) {
double tt, w;
double *b, *c; /* Base of parameters for this section */
int di = p->di;
int fdi = p->fdi;
int e, f, k;
double iparam = 0.0;
double oparam = 0.0;
double dd;
if (p->opt_msk & oc_i) {
dd = SHAPE_WEIGHT/(double)(di);
b = p->v + p->shp_off;
c = dav + p->shp_off;
for (e = 0; e < di; e++) {
for (k = 0; k < p->iluord[e]; k++) {
if (k <= 1) {
w = SHAPE_HW01;
} else if (k <= SHAPE_HBREAK) {
double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
w *= p->shp_smooth[e];
} else {
w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
w *= p->shp_smooth[e];
}
tt = *b++;
*c++ += 2.0 * dd * w * tt;
tt *= tt; /* Squared */
iparam += w * tt;
}
}
iparam *= dd;
}
if (p->opt_msk & oc_o) {
dd = SHAPE_WEIGHT/(double)(fdi);
b = p->v + p->out_off;
c = dav + p->out_off;
for (f = 0; f < fdi; f++) {
for (k = 0; k < p->oluord[f]; k++) {
if (k <= 1) {
w = SHAPE_HW01;
} else if (k <= SHAPE_HBREAK) {
double bl = (k - 1.0)/(SHAPE_HBREAK - 1.0);
w = (1.0 - bl) * SHAPE_HW01 + bl * SHAPE_HWBR;
w *= p->out_smooth[f];
} else {
w = SHAPE_HWBR + (k-SHAPE_HBREAK) * SHAPE_HWINC;
w *= p->out_smooth[f];
}
tt = *b++;
*c++ += 2.0 * dd * w * tt;
tt *= tt; /* Squared */
oparam += w * tt;
}
}
oparam *= dd;
}
return iparam + oparam;
}
/* Scale the shaper derivatives */
static void dshapscale(
xfit *p, /* Base of optimisation structure */
double *dav, /* del's */
double scale /* Scale factor */
) {
double tt, w;
double *b, *c; /* Base of parameters for this section */
int di = p->di;
int fdi = p->fdi;
int e, f, k;
if (p->opt_msk & oc_i) {
c = dav + p->shp_off;
for (e = 0; e < di; e++) {
for (k = 0; k < p->iluord[e]; k++) {
*c++ *= scale;
}
}
}
if (p->opt_msk & oc_o) {
c = dav + p->out_off;
for (f = 0; f < fdi; f++) {
for (k = 0; k < p->oluord[f]; k++) {
*c++ *= scale;
}
}
}
}
/* Progress function */
static void xfitprog(void *pdata, int perc) {
xfit *p = (xfit *)pdata;
if (p->verb) {
printf("%c% 3d%%",cr_char,perc);
if (perc == 100)
printf("\n");
fflush(stdout);
}
}
int xfitfunc_trace = 1;
/* Shaper+Matrix optimisation function handed to powell() */
/* We simply minimize the total delta E squared, consistent with smoothness */
static double xfitfunc(void *edata, double *v) {
xfit *p = (xfit *)edata;
double tw = 0.0; /* Total weight */
double ev = 0.0, rv, smv;
double tin[MXDI], out[MXDO];
int di = p->di;
int fdi = p->fdi;
int i, e, f;
/* Copy the parameters being optimised into xfit structure */
/* Special case - a single shaper curve. The first sm_iluord params */
/* are the common curve parameters, and the remainder are the matrix onwards */
if (p->opt_ssch) {
for (e = 0; e < di; e++) { /* Duplicate and extend to per channel curve params */
for (i = 0; i < p->sm_iluord; i++)
p->v[p->shp_offs[e] + i] = v[i];
for (; i < p->iluord[e]; i++)
p->v[p->shp_offs[e] + i] = 0.0;
}
for (i = p->sm_iluord; i < p->opt_cnt; i++)
p->v[p->mat_off + i - p->sm_iluord] = v[i];
} else {
for (i = 0; i < p->opt_cnt; i++) {
//printf("~1 param %d = %f\n",i,v[i]);
p->v[p->opt_off + i] = v[i];
}
}
/* For all our data points */
for (i = 0; i < p->nodp; i++) {
double del;
/* Apply input shaper channel curves */
for (e = 0; e < di; e++)
tin[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], p->rpoints[i].p[e],
p->in_min[e], p->in_max[e]);
/* Apply matrix cube interpolation */
icxCubeInterp(p->v + p->mat_off, fdi, di, out, tin);
/* Apply output channel curves */
for (f = 0; f < fdi; f++) {
if (p->flags & XFIT_OUT_LAB) {
out[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], out[f],
p->out_min[f], p->out_max[f]);
} else {
out[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], out[f],
p->out_min[f], p->out_max[f]);
}
}
/* Evaluate the error squared */
if (p->flags & XFIT_FM_INPUT) {
double pp[MXDI];
for (e = 0; e < di; e++)
pp[e] = p->rpoints[i].p[e];
for (f = 0; f < fdi; f++) {
double t1 = p->rpoints[i].v[f] - out[f]; /* Error in output */
/* Create input point offset by equivalent delta to output point */
/* error, in proportion to the partial derivatives for that output. */
for (e = 0; e < di; e++)
pp[e] += t1 * p->piv[i].ide[f][e];
}
del = p->to_de2(p->cntx2, pp, p->rpoints[i].p);
} else {
del = p->to_de2(p->cntx2, out, p->rpoints[i].v);
}
if (CURVEPOW > 1.0)
del = pow(del, CURVEPOW);
tw += p->rpoints[i].w;
ev += p->rpoints[i].w * del;
}
/* Normalise error to be an average delta E squared */
ev /= tw;
/* Sum with shaper parameters squared, to */
/* minimise unsconstrained "wiggles" */
smv = shapmag(p);
if (CURVEPOW > 1.0)
smv = pow(smv, CURVEPOW);
rv = ev + smv;
#ifdef DEBUG_PROGRESS
if (xfitfunc_trace)
fprintf(stdout,"~1(sm %f, ev %f)xfitfunc returning %f\n",smv,ev,rv);
#endif
return rv;
}
/* Shaper+Matrix optimisation function with partial derivatives, */
/* handed to conjgrad() */
static double dxfitfunc(void *edata, double *dv, double *v) {
xfit *p = (xfit *)edata;
double tw = 0.0; /* Total weight */
double ev = 0.0, rv, smv;
double tin[MXDI], out[MXDO];
double dav[MXPARMS]; /* Overall del due to del param vals */
double sdav[MXPARMS]; /* Overall del due to del smooth param vals */
double dtin_iv[MXDI * MXLUORD]; /* Del in itrans out due to del itrans param vals */
double dmato_mv[1 << MXDI]; /* Del in mat out due to del in matrix param vals */
double dmato_tin[MXDO * MXDI]; /* Del in mat out due to del in matrix input values */
double dout_ov[MXDO * MXLUORD]; /* Del in otrans out due to del in otrans param values */
double dout_mato[MXDO]; /* Del in otrans out due to del in otrans input values */
double dout_de[2][MXDIDO]; /* Del in DE due to two output values */
int di = p->di;
int fdi = p->fdi;
int i, jj, k, e, ee, f, ff;
/* Copy the parameters being optimised into xfit structure */
/* Special case - a single shaper curve. The first sm_iluord params */
/* are the common curve parameters, and the remainder are the matrix onwards */
if (p->opt_ssch) {
for (e = 0; e < di; e++) { /* Duplicate and extend to per channel curve params */
for (i = 0; i < p->sm_iluord; i++)
p->v[p->shp_offs[e] + i] = v[i];
for (; i < p->iluord[e]; i++)
p->v[p->shp_offs[e] + i] = 0.0;
}
for (i = p->sm_iluord; i < p->opt_cnt; i++)
p->v[p->mat_off + i - p->sm_iluord] = v[i];
} else {
for (i = 0; i < p->opt_cnt; i++) {
p->v[p->opt_off + i] = v[i];
}
}
/* Zero the accumulated partial derivatives */
/* We compute deriv for all parameters (not just current optimised) */
for (i = 0; i < p->tot_cnt; i++)
dav[i] = 0.0;
/* For all our data points */
for (i = 0; i < p->nodp; i++) {
double del;
/* Apply input channel curves */
for (e = 0; e < di; e++)
tin[e] = icxdpSTransFunc(p->v + p->shp_offs[e], &dtin_iv[p->shp_offs[e] - p->shp_off],
p->iluord[e], p->rpoints[i].p[e], p->in_min[e], p->in_max[e]);
/* Apply matrix cube interpolation */
icxdpdiCubeInterp(p->v + p->mat_off, dmato_mv, dmato_tin, fdi, di, out, tin);
/* Apply output channel curves */
for (f = 0; f < fdi; f++) {
if (p->flags & XFIT_OUT_LAB)
out[f] = icxdpdiSTransFunc(p->v + p->out_offs[f],
&dout_ov[p->out_offs[f] - p->out_off], &dout_mato[f],
p->oluord[f], out[f], p->out_min[f], p->out_max[f]);
else
out[f] = icxdpdiSTransFuncY2L(p->v + p->out_offs[f],
&dout_ov[p->out_offs[f] - p->out_off], &dout_mato[f],
p->oluord[f], out[f], p->out_min[f], p->out_max[f]);
}
/* Convert to Delta E and compute pde's into dout_de squared */
if (p->flags & XFIT_FM_INPUT) {
double tdout_de[2][MXDIDO];
double pp[MXDI];
for (e = 0; e < di; e++)
pp[e] = p->rpoints[i].p[e];
for (f = 0; f < fdi; f++) {
double t1 = p->rpoints[i].v[f] - out[f]; /* Error in output */
/* Create input point offset by equivalent delta to output point */
/* error, in proportion to the partial derivatives for that output. */
for (e = 0; e < di; e++)
pp[e] += t1 * p->piv[i].ide[f][e];
}
del = p->to_dde2(p->cntx2, tdout_de, pp, p->rpoints[i].p);
if (CURVEPOW > 1.0) {
double dadj;
dadj = CURVEPOW * pow(del, CURVEPOW - 1.0); /* Adjust derivative accordingly */
del = pow(del, CURVEPOW);
for (e = 0; e < di; e++)
tdout_de[0][e] *= dadj;
}
/* Compute partial derivative */
for (e = 0; e < di; e++) {
dout_de[0][e] = 0.0;
for (f = 0; f < fdi; f++) {
dout_de[0][e] += tdout_de[0][e] * p->piv[i].ide[f][e];
}
}
} else {
del = p->to_dde2(p->cntx2, dout_de, out, p->rpoints[i].v);
if (CURVEPOW > 1.0) {
double dadj;
dadj = CURVEPOW * pow(del, CURVEPOW - 1.0); /* Adjust derivative accordingly */
del = pow(del, CURVEPOW);
for (f = 0; f < fdi; f++)
dout_de[0][f] *= dadj;
}
}
/* Accumulate total weighted delta E squared */
tw += p->rpoints[i].w;
ev += p->rpoints[i].w * del;
/* Compute and accumulate partial difference values for each parameter value */
if (p->opt_msk & oc_i) {
/* Input transfer parameters */
for (ee = 0; ee < di; ee++) { /* Parameter input chanel */
for (k = 0; k < p->iluord[ee]; k++) { /* Param within channel */
double vv = 0.0;
jj = p->shp_offs[ee] - p->shp_off + k; /* Overall input trans param */
// for (ff = 0; ff < 3; ff++) { /* Lab channels */
for (ff = 0; ff < fdi; ff++) { /* Output channels */
vv += dout_de[0][ff] * dout_mato[ff]
* dmato_tin[ff * di + ee] * dtin_iv[jj];
}
dav[p->shp_off + jj] += p->rpoints[i].w * vv;
}
}
}
if (p->opt_msk & oc_m) {
/* Matrix parameters */
for (ff = 0; ff < fdi; ff++) { /* Parameter output chanel */
for (ee = 0; ee < (1 << di); ee++) { /* Matrix input combination chanel */
double vv = 0.0;
jj = ff * (1 << di) + ee; /* Matrix Parameter index */
vv += dout_de[0][ff] * dout_mato[ff] * dmato_mv[ee];
dav[p->mat_off + jj] += p->rpoints[i].w * vv;
}
}
}
if (p->opt_msk & oc_o) {
/* Output transfer parameters */
for (ff = 0; ff < fdi; ff++) { /* Parameter output chanel */
for (k = 0; k < p->oluord[ff]; k++) { /* Param within channel */
double vv = 0.0;
jj = p->out_offs[ff] - p->out_off + k; /* Overall output trans param */
vv += dout_de[0][ff] * dout_ov[jj];
dav[p->out_off + jj] += p->rpoints[i].w * vv;
}
}
}
}
/* Normalise error to be an average delta E squared */
ev /= tw;
for (i = 0; i < p->tot_cnt; i++) {
dav[i] /= tw;
sdav[i] = 0.0;
}
/* Sum with shaper parameters squared, to */
/* minimise unsconstrained "wiggles" */
/* Compute partial derivative wrt those parameters too */
smv = dshapmag(p, sdav);
if (CURVEPOW > 1.0) {
double dadj;
dadj = CURVEPOW * pow(smv, CURVEPOW - 1.0); /* Adjust derivative accordingly */
smv = pow(smv, CURVEPOW);
dshapscale(p, sdav, dadj); /* Scale the partial derivatives */
}
rv = ev + smv;
/* Sum the del for parameters being optimised and copy to return array */
if (p->opt_ssch) {
for (i = 0; i < p->sm_iluord; i++)
dv[i] = 0.0;
for (e = 0; e < di; e++) { /* Combine per channel curve de's */
for (i = 0; i < p->sm_iluord; i++)
dv[i] += dav[p->shp_offs[e] + i] + sdav[p->shp_offs[e] + i];
}
for (i = p->sm_iluord; i < p->opt_cnt; i++) /* matrix and rest de's */
dv[i] = dav[p->mat_off + i - p->sm_iluord] + sdav[p->mat_off + i - p->sm_iluord];
} else {
for (i = 0; i < p->opt_cnt; i++)
dv[i] = dav[p->opt_off + i] + sdav[p->opt_off + i];
}
#ifdef DEBUG_PROGRESS
fprintf(stdout,"~1(sm %f, ev %f)dxfitfunc returning %f\n",smv,ev,rv);
#endif
return rv;
}
#ifdef NEVER
/* Check partial derivative function within xfitfunc() [Intensive check] */
static double _xfitfunc(void *edata, double *v) {
xfit *p = (xfit *)edata;
int i;
double dv[MXPARMS];
double rv, drv;
double trv;
int verb;
rv = xfitfunc(edata, v);
verb = p->verb;
p->verb = 0;
drv = dxfitfunc(edata, dv, v);
p->verb = verb;
if (fabs(rv - drv) > 1e-6)
printf("######## RV MISMATCH is %f should be %f ########\n",rv,drv);
/* Check each parameter delta */
xfitfunc_trace = 0;
for (i = 0; i < p->opt_cnt; i++) {
double del;
v[i] += 1e-7;
trv = xfitfunc(edata, v);
v[i] -= 1e-7;
/* Check that del is correct */
del = (trv - rv)/1e-7;
if (fabs(dv[i] - del) > 0.04) {
//printf("~1 del = %f from (trv %f - rv %f)/0.1\n",del,trv,rv);
printf("######## EXCESSIVE at v[%d] is %f should be %f ########\n",i,dv[i],del);
}
}
xfitfunc_trace = 1;
return rv;
}
#define xfitfunc _xfitfunc
#endif /* NEVER */
#ifdef NEVER
/* Check partial derivative function within dxfitfunc() [Less intensive check] */
static double _dxfitfunc(void *edata, double *dv, double *v) {
xfit *p = (xfit *)edata;
int i;
double rv, drv;
double trv;
int verb;
int exec = 0;
rv = xfitfunc(edata, v);
verb = p->verb;
p->verb = 0;
drv = dxfitfunc(edata, dv, v);
p->verb = verb;
if (fabs(rv - drv) > 1e-6)
printf("######## RV MISMATCH is %f should be %f ########\n",rv,drv);
/* Check each parameter delta */
xfitfunc_trace = 0;
for (i = 0; i < p->opt_cnt; i++) {
double del;
v[i] += 1e-7;
trv = xfitfunc(edata, v);
v[i] -= 1e-7;
/* Check that del is correct */
del = (trv - rv)/1e-7;
if (fabs(dv[i] - del) > 0.04) {
//printf("~1 del = %f from (trv %f - rv %f)/0.1\n",del,trv,rv);
printf("######## EXCESSIVE at v[%d] is %f should be %f ########\n",i,dv[i],del);
exec = 1;
}
}
#ifdef NEVER
if (exec) {
printf("~1 parameters are:\n");
for (i = 0; i < p->opt_cnt; i++)
printf("p->wv[%d] = %f;\n",i,v[i]);
exit(1);
}
#endif
xfitfunc_trace = 1;
return rv;
}
#define dxfitfunc _dxfitfunc
#endif /* NEVER */
/* - - - - - - - - - */
/* Output curve symetry optimisation function handed to powell() */
/* Just the order 0 value will be adjusted */
static double symoptfunc(void *edata, double *v) {
xfit *p = (xfit *)edata;
double out[1], in[1] = { 0.0 };
int ch = p->opt_ch; /* Output channel being adjusted for symetry */
double rv;
/* Copy the parameter being tested back into xfit */
p->v[p->out_offs[ch]] = v[0];
if (p->flags & XFIT_OUT_LAB)
*out = icxSTransFunc(p->v + p->out_offs[ch], p->oluord[ch], *in,
p->out_min[ch], p->out_max[ch]);
else
*out = icxSTransFuncY2L(p->v + p->out_offs[ch], p->oluord[ch], *in,
p->out_min[ch], p->out_max[ch]);
rv = out[0] * out[0];
#ifdef DEBUG_PROGRESS
printf("~1symoptfunc returning %f\n",rv);
#endif
return rv;
}
/* - - - - - - - - - */
/* Set up for an optimisation run: */
/* Figure out parameters being optimised, */
/* copy them to start values, */
/* init and scale the search radius */
static void setup_xfit(
xfit *p,
double *wv, /* Return parameters to hand to optimiser */
double *sa, /* Return search radius to hand to optimiser */
double transrad,/* Nominal transfer curve radius, 0.0 - 3.0 */
double matrad /* Nominal matrix radius, 0.0 - 1.0 */
) {
int i;
p->opt_off = -1;
p->opt_cnt = 0;
if (p->opt_msk & oc_i) {
if (p->opt_ssch) { /* Special case - should only be used first, */
/* Fitting a sigle common input shaper curve. */
if (p->opt_off < 0)
p->opt_off = p->mat_off - p->sm_iluord; /* Shouldn't be used... */
p->opt_cnt += p->sm_iluord;
for (i = 0; i < p->sm_iluord; i++) {
*wv++ = 0.0;
*sa++ = transrad;
}
} else { /* Initial or continuing fitting of all the curves */
if (p->opt_off < 0)
p->opt_off = p->shp_off;
p->opt_cnt += p->shp_cnt;
for (i = 0; i < p->shp_cnt; i++) {
*wv++ = p->v[p->shp_off + i];
*sa++ = transrad;
}
}
}
if (p->opt_msk & oc_m) {
if (p->opt_off < 0)
p->opt_off = p->mat_off;
p->opt_cnt += p->mat_cnt;
for (i = 0; i < p->mat_cnt; i++) {
*wv++ = p->v[p->mat_off + i];
*sa++ = matrad;
}
}
if (p->opt_msk & oc_o) {
if (p->opt_off < 0)
p->opt_off = p->out_off;
p->opt_cnt += p->out_cnt;
for (i = 0; i < p->out_cnt; i++) {
*wv++ = p->v[p->out_off + i];
*sa++ = transrad;
}
}
if (p->opt_cnt > MXPARMS)
error("setup_xfit: asert, %d exceeded MXPARMS %d",p->opt_cnt,MXPARMS);
#ifdef DEBUG
printf("setup_xfit() got opt_msk 0x%x, opt_off %d, opt_cnt = %d\n",
p->opt_msk,p->opt_off,p->opt_cnt);
#endif /* DEBUG */
}
#ifdef DEBUG
/* Diagnostic */
static void dump_xfit(
xfit *p
) {
int i, e, f;
double *b; /* Base of parameters for this section */
int di, fdi;
di = p->di;
fdi = p->fdi;
/* Input positioning curve */
b = p->v + p->pos_off;
for (e = 0; e < di; b += p->iluord[e], e++) {
printf("pos %d = ",e);
for (i = 0; i < p->iluord[e]; i++)
printf("%f ",b[i]);
printf("\n");
}
/* Input shaper curve */
b = p->v + p->shp_off;
for (e = 0; e < di; b += p->iluord[e], e++) {
printf("shp %d = ",e);
for (i = 0; i < p->iluord[e]; i++)
printf("%f ",b[i]);
printf("\n");
}
/* Matrix */
b = p->v + p->mat_off;
for (e = 0; e < (1 << di); e++) {
printf("mx %d = ",e);
for (f = 0; f < fdi; f++)
printf("%f ",*b++);
printf("\n");
}
/* Output curve */
b = p->v + p->out_off;
for (f = 0; f < fdi; b += p->oluord[f], f++) {
printf("out %d = ",f);
for (i = 0; i < p->oluord[f]; i++)
printf("%f ",b[i]);
printf("\n");
}
}
#endif /* DEBUG */
/* - - - - - - - - - */
/* Setup the pseudo inverse information for each test point, */
/* using the current model. */
static void setup_piv(xfit *p) {
int di = p->di;
int fdi = p->fdi;
int i, e, f;
/* Estimate in -> out partial derivatives */
for (i = 0; i < p->nodp; i++) {
double pd[MXDO][MXDI];
double pp[MXDI];
double vv[MXDIDO];
/* Estimate in -> out partial derivatives */
for (e = 0; e < di; e++)
pp[e] = p->ipoints[i].p[e];
/* Apply input shaper channel curves */
for (e = 0; e < di; e++)
vv[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], pp[e],
p->in_min[e], p->in_max[e]);
/* Apply matrix cube interpolation */
icxCubeInterp(p->v + p->mat_off, fdi, di, vv, vv);
/* Apply output channel curves */
for (f = 0; f < fdi; f++) {
if (p->flags & XFIT_OUT_LAB)
vv[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], vv[f],
p->out_min[f], p->out_max[f]);
else
vv[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], vv[f],
p->out_min[f], p->out_max[f]);
}
for (e = 0; e < di; e++) {
double tt[MXDIDO];
pp[e] += 1e-4;
/* Apply input shaper channel curves */
for (e = 0; e < di; e++)
tt[e] = icxSTransFunc(p->v + p->shp_offs[e], p->iluord[e], pp[e],
p->in_min[e], p->in_max[e]);
/* Apply matrix cube interpolation */
icxCubeInterp(p->v + p->mat_off, fdi, di, tt, tt);
/* Apply output channel curves */
for (f = 0; f < fdi; f++) {
if (p->flags & XFIT_OUT_LAB)
tt[f] = icxSTransFunc(p->v + p->out_offs[f], p->oluord[f], tt[f],
p->out_min[f], p->out_max[f]);
else
tt[f] = icxSTransFuncY2L(p->v + p->out_offs[f], p->oluord[f], tt[f],
p->out_min[f], p->out_max[f]);
}
for (f = 0; f < p->fdi; f++)
pd[f][e] = (tt[f] - vv[f])/1e-4;
pp[e] -= 1e-4;
}
/* Compute a psudo inverse matrix to map rout delta E to */
/* in delta E in proportion to the pd magnitude. */
for (f = 0; f < fdi; f++) {
double ss = 0.0;
for (e = 0; e < di; e++) /* Sum of pd's ^4 */
ss += pd[f][e] * pd[f][e] * pd[f][e] * pd[f][e];
ss = sqrt(ss);
if (ss > 1e-8) {
for (e = 0; e < di; e++)
p->piv[i].ide[f][e] = pd[f][e]/ss;
} else { /* Hmm. */
for (e = 0; e < di; e++)
p->piv[i].ide[f][e] = 0.0;
}
}
}
}
/* - - - - - - - - - */
/* Function to pass to rspl to re-set output values, */
/* to account for skeleton model offset. */
static void
skm_rspl_out(
void *pp, /* relativectx structure */
double *out, /* output value */
double *in /* input value */
) {
xfit *p = (xfit *)pp;
int f, fdi = p->fdi;
double inval[MXDI];
double skval[MXDO];
/* Look up the skeleton value for this grid point */
xfit_invinpscurves(p, inval, in); /* Back to input values */
p->skm->lookup(p->skm, skval, inval); /* Skm */
xfit_abs_to_rel(p, skval, skval);
xfit_invoutcurves(p, skval, skval);
for (f = 0; f < fdi; f++)
out[f] += skval[f]; /* Add it back */
}
/* Weak function rspl callback (not used) */
void skm_weak(void *cbntx, double *out, double *in) {
xfit *p = (xfit *)cbntx;
#ifndef NEVER
int f, fdi = p->fdi;
for (f = 0; f < fdi; f++)
out[f] = 0.0; /* Deviation from skeleton should tend to zero */
#else /* Skeleton as weak atractor */
int f, fdi = p->fdi;
double inval[MXDI];
/* Look up the skeleton value for this grid point */
xfit_invinpscurves(p, inval, in); /* Back to input values */
p->skm->lookup(p->skm, out, inval); /* Skm */
xfit_abs_to_rel(p, out, out);
xfit_invoutcurves(p, out, out);
#endif
}
/* - - - - - - - - - */
/* Function to pass to rspl to re-set output values, */
/* to make them relative to the white and black points */
static void
conv_rspl_out(
void *pp, /* relativectx structure */
double *out, /* output value */
double *in /* input value */
) {
xfit *p = (xfit *)pp;
double tt[3];
/* Convert the clut values to output values */
xfit_outcurves(p, tt, out);
if (p->flags & XFIT_OUT_LAB) {
icmLab2XYZ(&icmD50, tt, tt);
icmMulBy3x3(out, p->cmat, tt);
icmXYZ2Lab(&icmD50, out, out);
} else { /* We are all in XYZ */
icmMulBy3x3(out, p->cmat, tt);
}
/* And then convert them back to clut values */
xfit_invoutcurves(p, out, out);
}
/* Function to pass to rspl to re-set output values, */
/* to clip any with Y over 1.0 to D50 */
static void
clip_rspl_out(
void *pp, /* relativectx structure */
double *out, /* output value */
double *in /* input value */
) {
xfit *p = (xfit *)pp;
double tt[3];
/* Convert the clut values to output values */
xfit_outcurves(p, tt, out);
if (p->flags & XFIT_OUT_LAB) {
if (tt[0] > 100.0)
icmCpy3(out, p->cmat[0]);
} else {
if (tt[1] > 1.0)
icmCpy3(out, p->cmat[0]);
}
}
//#ifdef SPECIAL_TEST
/* - - - - - - - - - */
/* Execute the linear XYZ device model */
static void domodel(double *out, double *in) {
double tmp[3];
int i, j;
double col[3][3]; /* sRGB additive colorant values in XYZ :- [out][in] */
col[0][0] = 0.412424; /* X from R */
col[0][1] = 0.357579; /* X from G */
col[0][2] = 0.180464; /* X from B */
col[1][0] = 0.212656; /* Y from R */
col[1][1] = 0.715158; /* Y from G */
col[1][2] = 0.0721856; /* Y from B */
col[2][0] = 0.0193324; /* Z from R */
col[2][1] = 0.119193; /* Z from G */
col[2][2] = 0.950444; /* Z from B */
#ifdef SPECIAL_TEST_GAMMA
tmp[0] = pow(in[0], 1.9);
tmp[1] = pow(in[1], 2.0);
tmp[2] = pow(in[2], 2.1);
#else
tmp[0] = in[0];
tmp[1] = in[1];
tmp[2] = in[2];
#endif
for (j = 0; j < 3; j++) {
out[j] = 0.0;
for (i = 0; i < 3; i++)
out[j] += col[j][i] * tmp[i];
}
}
#ifdef SPECIAL_FORCE
/* Function to pass to rspl to set nodes against */
/* synthetic model. */
static void
set_rspl_out1(
void *pp, /* relativectx structure */
double *out, /* output value */
double *in /* input value */
) {
xfit *p = (xfit *)pp;
double tt[3], tout[3];
/* Convert the input' values to input values */
xfit_invinpscurves(p, tt, in);
/* Synthetic linear rgb->XYZ model */
domodel(tout, tt);
/* Apply abs->rel white point adjustment */
icmMulBy3x3(tout, p->cmat, tout);
#ifdef SPECIAL_TEST_LAB
icmXYZ2Lab(&icmD50, tout, tout);
#endif
/* And then convert them back to clut values */
xfit_invoutcurves(p, tout, tout);
#ifdef DEBUG
printf("~1 changing %f %f %f -> %f %f %f\n", out[0], out[1], out[2], tout[0], tout[1], tout[2]);
#endif
out[0] = tout[0];
out[1] = tout[1];
out[2] = tout[2];
}
#endif /* SPECIAL_FORCE */
/* - - - - - - - - - */
/* Do the fitting. */
/* return nz on error */
/* 1 = malloc or other error */
static int xfit_fit(
struct _xfit *p,
int flags, /* Flag values */
int di, /* Input dimensions */
int fdi, /* Output dimensions */
int rsplflags, /* clut rspl creation flags */
double *wp, /* if flags & XFIT_OUT_WP_REL or XFIT_OUT_WP_REL_US, */
/* Initial white point, returns final wp */
double *dw, /* Device white value to adjust to be D50 */
double wpscale, /* If >= 0.0 scale final wp */
double *dgw, /* Device space gamut boundary white for XFIT_OUT_WP_REL_US */
/* (ie. RGB 1,1,1 CMYK 0,0,0,0, etc) */
cow *ipoints, /* Array of data points to fit - referece taken */
int nodp, /* Number of data points */
icxMatrixModel *skm, /* Optional skeleton model (used for input profiles) */
double in_min[MXDI], /* Input value scaling/domain minimum */
double in_max[MXDI], /* Input value scaling/domain maximum */
int gres[MXDI], /* clut resolutions being optimised for/returned */
double out_min[MXDO], /* Output value scaling/range minimum */
double out_max[MXDO], /* Output value scaling/range maximum */
// co *bpo, /* If != NULL, black point override in same spaces as ipoints */
double smooth, /* clut rspl smoothing factor */
double oavgdev[MXDO], /* Average output value deviation */
double demph, /* dark emphasis factor for cLUT grid res. */
int iord[], /* Order of input pos/shaper curve for each dimension */
int sord[], /* Order of input sub-grid shaper curve (not used) */
int oord[], /* Order of output shaper curve for each dimension */
double shp_smooth[MXDI],/* Smoothing factors for each curve, nom = 1.0 */
double out_smooth[MXDO],
optcomb tcomb, /* Flag - target elements to fit. */
void *cntx2, /* Context of callbacks */
/* Callback to convert two fit values delta E squared */
double (*to_de2)(void *cntx, double *in1, double *in2),
/* Same as above, with partial derivatives */
double (*to_dde2)(void *cntx, double dout[2][MXDIDO], double *in1, double *in2)
) {
int i, e, f;
double *b; /* Base of parameters for this section */
int poff;
double powtol = POWTOL1; /* powell/conjgrad initial tollerance */
int maxits = MAXITS1; /* powell/conjgrad initial maximum itterations */
if (tcomb & oc_io) /* If we're doing anything, we need the matrix */
tcomb |= oc_m;
p->flags = flags;
if (flags & XFIT_VERB)
p->verb = 1;
else
p->verb = 0;
p->di = di;
p->fdi = fdi;
p->wp = wp; /* Take reference, so modified wp can be returned */
p->dw = dw;
p->nodp = nodp;
p->skm = skm; /* This isn't current used by profin, because it doesn't help.. */
p->ipoints = ipoints;
p->tcomb = tcomb;
p->cntx2 = cntx2;
for (e = 0; e < di; e++)
p->gres[e] = gres[e];
p->to_de2 = to_de2;
p->to_dde2 = to_dde2;
#ifdef DEBUG
printf("xfit_fit called with flags = 0x%x, di = %d, fdi = %d, nodp = %d, tcomb = 0x%x\n",flags,di,fdi,nodp,tcomb);
#endif
//printf("~1 out min = %f %f %f max = %f %f %f\n", out_min[0], out_min[1], out_min[2], out_max[0], out_max[1], out_max[2]);
/* Sanity protect shaper orders */
/* and save scaling and smoothness factors. */
p->sm_iluord = MXLUORD+1;
for (e = 0; e < di; e++) {
if (iord[e] > MXLUORD)
p->iluord[e] = MXLUORD;
else
p->iluord[e] = iord[e];
if (p->iluord[e] < p->sm_iluord)
p->sm_iluord = p->iluord[e];
p->in_min[e] = in_min[e];
p->in_max[e] = in_max[e];
p->shp_smooth[e] = shp_smooth[e];
}
for (f = 0; f < fdi; f++) {
if (oord[f] > MXLUORD)
p->oluord[f] = MXLUORD;
else
p->oluord[f] = oord[f];
p->out_min[f] = out_min[f];
p->out_max[f] = out_max[f];
p->out_smooth[f] = out_smooth[f];
}
/* Compute parameter offset and count information */
p->shp_off = 0;
for (poff = p->shp_off, p->shp_cnt = 0, e = 0; e < di; e++) {
p->shp_offs[e] = poff;
p->shp_cnt += p->iluord[e];
poff += p->iluord[e];
}
p->mat_off = p->shp_off + p->shp_cnt;
for (poff = p->mat_off, p->mat_cnt = 0, f = 0; f < fdi; f++) {
p->mat_offs[f] = poff;
p->mat_cnt += (1 << di);
poff += (1 << di);
}
p->out_off = p->mat_off + p->mat_cnt;
for (poff = p->out_off, p->out_cnt = 0, f = 0; f < fdi; f++) {
p->out_offs[f] = poff;
p->out_cnt += p->oluord[f];
poff += p->oluord[f];
}
p->pos_off = p->out_off + p->out_cnt;
for (poff = p->pos_off, p->pos_cnt = 0, e = 0; e < di; e++) {
p->pos_offs[e] = poff;
p->pos_cnt += p->iluord[e];
poff += p->iluord[e];
}
p->tot_cnt = p->shp_cnt + p->mat_cnt + p->out_cnt + p->pos_cnt;
if (p->tot_cnt > MXPARMS)
error("xfit_fit: assert tot_cnt exceeds MXPARMS");
/* Allocate space for parameter values */
if (p->v != NULL) {
free(p->v);
p->v = NULL;
}
if (p->wv != NULL) {
free(p->wv);
p->wv = NULL;
}
if (p->sa != NULL) {
free(p->sa);
p->sa = NULL;
}
if ((p->v = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
return 1;
if ((p->wv = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
return 1;
if ((p->sa = (double *)calloc(p->tot_cnt, sizeof(double))) == NULL)
return 1;
/* Setup initial white point abs->rel conversions */
if ((p->flags & XFIT_OUT_WP_REL) != 0) {
icmXYZNumber _wp;
icmAry2XYZ(_wp, p->wp);
/* Absolute->Aprox. Relative Adaptation matrix and */
/* Aproximate relative to absolute conversion matrix */
if (p->picc != NULL) {
p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, p->toAbs, p->fromAbs, icmD50, _wp);
} else {
icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
}
if (p->verb) {
double lab[3];
icmXYZ2Lab(&icmD50, lab, p->wp);
printf("Initial White Point XYZ %f %f %f, Lab %f %f %f\n",
p->wp[0], p->wp[1], p->wp[2], lab[0], lab[1], lab[2]);
}
} else {
icmSetUnity3x3(p->fromAbs);
icmSetUnity3x3(p->toAbs);
icmCpy3(p->wp, icmD50_ary3); /* Default to forced absolute */
}
/* Setup input position/shape curves to be linear initially */
b = p->v + p->shp_off;
for (e = 0; e < di; b += p->iluord[e], e++) {
for (i = 0; i < p->iluord[e]; i++) {
b[i] = 0.0;
}
}
/* Setup matrix to be pure colorant' values initially */
b = p->v + p->mat_off;
for (e = 0; e < (1 << di); e++) { /* For each colorant combination */
int j, k, bk = 0;
double bdif = 1e6;
double ov[MXDO];
/* Search the patch list to find the one closest to this input combination */
for (k = 0; k < p->nodp; k++) {
double dif = 0.0;
for (j = 0; j < di; j++) {
double tt;
if (e & (1 << j))
tt = p->in_max[j] - p->ipoints[k].p[j];
else
tt = p->in_min[j] - p->ipoints[k].p[j];
dif += tt * tt;
}
if (dif < bdif) { /* best so far */
bdif = dif;
bk = k;
if (dif < 0.001)
break; /* Don't bother looking further */
}
}
xfit_abs_to_rel(p, ov, p->ipoints[bk].v);
for (f = 0; f < fdi; f++)
b[f * (1 << di) + e] = ov[f];
}
/* Setup output curves to be linear initially */
b = p->v + p->out_off;
for (f = 0; f < fdi; b += p->oluord[f], f++) {
for (i = 0; i < p->oluord[f]; i++) {
b[i] = 0.0;
}
}
/* Setup positioning curves to be linear initially */
b = p->v + p->pos_off;
for (e = 0; e < di; b += p->iluord[e], e++) {
for (i = 0; i < p->iluord[e]; i++) {
b[i] = 0.0;
}
}
/* Create copy of input points with output converted to white relative */
if (p->rpoints == NULL) {
if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
return 1;
}
for (i = 0; i < p->nodp; i++) {
p->rpoints[i].w = p->ipoints[i].w;
for (e = 0; e < di; e++)
p->rpoints[i].p[e] = p->ipoints[i].p[e];
for (f = 0; f < fdi; f++)
p->rpoints[i].v[f] = p->ipoints[i].v[f];
/* out -> rout */
xfit_abs_to_rel(p, p->rpoints[i].v, p->rpoints[i].v);
}
/* Allocate array of pseudo-inverse matricies */
if ((p->flags & XFIT_FM_INPUT) != 0 && p->piv == NULL) {
if ((p->piv = (xfit_piv *)malloc(p->nodp * sizeof(xfit_piv))) == NULL)
return 1;
}
/* Allocate array of span DE's for current opt channel */
{
int lres = 0;
for (e = 0; e < di; e++) {
if (p->gres[e] > lres)
lres = p->gres[e];
}
if ((p->uerrv = (double *)malloc(lres * sizeof(double))) == NULL)
return 1;
}
/* Do the fitting one part at a time, then together */
/* Shaper curves are created if position or shaper curves are requested */
/* Fit just the matrix */
if ((p->tcomb & oc_ipo) != 0
&& (p->tcomb & oc_m) == oc_m) { /* Only bother with matrix if in and/or out */
double rerr;
if (p->verb)
printf("About to optimise temporary matrix\n");
/* Setup pseudo-inverse if we need it */
if (p->flags & XFIT_FM_INPUT)
setup_piv(p);
#ifdef DEBUG
printf("\nBefore matrix opt:\n");
dump_xfit(p);
#endif
/* Optimise matrix on its own */
p->opt_ssch = 0;
p->opt_ch = -1;
p->opt_msk = oc_m;
setup_xfit(p, p->wv, p->sa, 0.0, 0.5);
#ifdef NODDV
if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#else
if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
warning("xfit_fit: Conjgrad failed to converge, residual error = %f", rerr);
#endif
for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
p->v[p->opt_off + i] = p->wv[i];
#ifdef DEBUG
printf("\nAfter matrix opt:\n");
dump_xfit(p);
#endif
}
/* Optimise input and matrix together */
if ((p->tcomb & oc_im) == oc_im) {
double rerr;
int sm_iluord = p->sm_iluord;
if (p->verb)
printf("About to optimise a common ord 0 input curve and matrix\n");
/* Setup pseudo-inverse if we need it */
if (p->flags & XFIT_FM_INPUT)
setup_piv(p);
p->opt_ssch = 1;
p->sm_iluord = 1; /* Do a single order for first up */
p->opt_ch = -1;
p->opt_msk = oc_im;
setup_xfit(p, p->wv, p->sa, 0.5, 0.3);
#ifdef NODDV
if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#endif
}
#else
if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
}
#endif /* !NODDV */
for (e = 0; e < di; e++) { /* Copy optimised values back */
for (i = 0; i < p->sm_iluord; i++)
p->v[p->shp_offs[e] + i] = p->wv[i];
for (; i < p->iluord[e]; i++)
p->v[p->shp_offs[e] + i] = 0.0;
}
for (i = p->sm_iluord; i < p->opt_cnt; i++)
p->v[p->mat_off + i - p->sm_iluord] = p->wv[i];
#ifdef DEBUG
printf("\nAfter single input and matrix opt:\n");
dump_xfit(p);
#endif
/* - - - - - - - - - - - */
if (p->verb)
printf("About to optimise a common input curve and matrix\n");
/* Setup pseudo-inverse if we need it */
if (p->flags & XFIT_FM_INPUT)
setup_piv(p);
p->opt_ssch = 1;
p->sm_iluord = sm_iluord; /* restore this */
p->opt_ch = -1;
p->opt_msk = oc_im;
setup_xfit(p, p->wv, p->sa, 0.5, 0.3);
#ifdef NODDV
if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#endif
}
#else
if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
}
#endif /* !NODDV */
for (e = 0; e < di; e++) { /* Copy optimised values back */
for (i = 0; i < p->sm_iluord; i++)
p->v[p->shp_offs[e] + i] = p->wv[i];
for (; i < p->iluord[e]; i++)
p->v[p->shp_offs[e] + i] = 0.0;
}
for (i = p->sm_iluord; i < p->opt_cnt; i++)
p->v[p->mat_off + i - p->sm_iluord] = p->wv[i];
#ifdef DEBUG
printf("\nAfter single input and matrix opt:\n");
dump_xfit(p);
#endif
/* - - - - - - - - - - - */
if (p->verb)
printf("About to optimise input curves and matrix\n");
if ((p->tcomb & oc_mo) != oc_mo) { /* If this will be last fit */
powtol = POWTOL;
maxits = MAXITS;
}
/* Setup pseudo-inverse if we need it */
if (p->flags & XFIT_FM_INPUT)
setup_piv(p);
p->opt_ssch = 0;
p->opt_ch = -1;
p->opt_msk = oc_im;
setup_xfit(p, p->wv, p->sa, 0.5, 0.3);
/* Suppress the warnings the first time through - it's better to cut off the */
/* itterations and move on to the output curve, and worry about it not */
/* converging the second time through. */
#ifdef NODDV
if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#endif
}
#else
if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0) {
#ifdef DEBUG
warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
}
#endif /* !NODDV */
for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
p->v[p->opt_off + i] = p->wv[i];
#ifdef DEBUG
printf("\nAfter input and matrix opt:\n");
dump_xfit(p);
#endif
}
/* Optimise the matrix and output curves together */
if ((p->tcomb & oc_mo) == oc_mo) {
double rerr;
if (p->verb)
printf("About to optimise output curves and matrix\n");
if ((p->tcomb & oc_im) != oc_im) { /* If this will be last fit */
powtol = POWTOL;
maxits = MAXITS;
}
/* Setup pseudo-inverse if we need it */
if (p->flags & XFIT_FM_INPUT)
setup_piv(p);
p->opt_ssch = 0;
p->opt_ch = -1;
p->opt_msk = oc_mo;
setup_xfit(p, p->wv, p->sa, 0.3, 0.3);
#ifdef NODDV
if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#else
if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits, xfitfunc,
dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
p->v[p->opt_off + i] = p->wv[i];
#ifdef DEBUG
printf("\nAfter output opt:\n");
dump_xfit(p);
#endif
/* Optimise input and matrix together again, after altering matrix */
if ((p->tcomb & oc_im) == oc_im) {
if (p->verb)
printf("About to optimise input curves and matrix again\n");
#ifndef NODDV
if ((p->tcomb & oc_imo) != oc_imo) /* If this will be last fit */
#endif
{
powtol = POWTOL;
maxits = MAXITS;
}
/* Setup pseudo-inverse if we need it */
if (p->flags & XFIT_FM_INPUT)
setup_piv(p);
p->opt_ssch = 0;
p->opt_ch = -1;
p->opt_msk = oc_im;
setup_xfit(p, p->wv, p->sa, 0.2, 0.2);
#ifdef NODDV
if (powell(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, (void *)p, xfitprog, (void *)p) != 0)
warning("xfit_fit: Powell failed to converge, residual error = %f",rerr);
#else
if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, powtol, maxits,
xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
#endif
for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
p->v[p->opt_off + i] = p->wv[i];
}
#ifdef DEBUG
printf("\nAfter 2nd input and matrix opt:\n");
dump_xfit(p);
#endif
#ifndef NODDV
/* Optimise all together */
/* (This is very slow using powell) */
if ((p->tcomb & oc_imo) == oc_imo) {
if (p->verb)
printf("About to optimise input, matrix and output together\n");
/* Setup pseudo-inverse if we need it */
if (p->flags & XFIT_FM_INPUT)
setup_piv(p);
p->opt_ssch = 0;
p->opt_ch = -1;
p->opt_msk = oc_imo;
setup_xfit(p, p->wv, p->sa, 0.1, 0.1);
if (conjgrad(&rerr, p->opt_cnt, p->wv, p->sa, POWTOL, MAXITS,
xfitfunc, dxfitfunc, (void *)p, xfitprog, (void *)p) != 0)
warning("xfit_fit: Conjgrad failed to converge, residual error = %f",rerr);
for (i = 0; i < p->opt_cnt; i++) /* Copy optimised values back */
p->v[p->opt_off + i] = p->wv[i];
/* Setup final pseudo-inverse from shaper/matrix/out */
if (p->flags & XFIT_FM_INPUT)
setup_piv(p);
}
#ifdef DEBUG
printf("\nAfter all together opt:\n");
dump_xfit(p);
#endif
#endif /* !NODDV */
/* Adjust output curve white point. */
/* This is for the benefit of the B2A table */
if (p->flags & XFIT_OUT_ZERO) {
if (p->verb)
printf("About to adjust a and b output curves for white point\n");
for (f = 1; f < 3 && f < p->fdi; f++) {
p->opt_ch = f;
p->wv[0] = p->v[p->out_offs[f]]; /* Current parameter value */
p->sa[0] = 0.1; /* Search radius */
if (powell(&rerr, 1, p->wv, p->sa, 0.0000001, 1000,
symoptfunc, (void *)p, NULL, NULL) != 0)
error("xfit_fit: Powell failed to converge, residual error = %f",rerr);
p->v[p->out_offs[f]] = p->wv[0]; /* Copy results back */
}
}
}
/* In case we don't generate position curves, */
/* copy the input curves to the position, so that */
/* ipos is computed correctly */
{
double *bb;
b = p->v + p->shp_off;
bb = p->v + p->pos_off;
for (e = 0; e < di; b += p->iluord[e], bb += p->iluord[e], e++) {
for (i = 0; i < p->iluord[e]; i++) {
bb[i] = b[i];
}
}
}
/* If we want position curves, generate them */
/* (This could possibly be improved by using some sort */
/* of optimization drivel approach rather than the predictive */
/* method used here.) */
if (p->tcomb & oc_p) {
int ee;
if (p->verb)
printf("About to create grid position input curves\n");
/* Allocate in->rout duplicate point set */
if (p->rpoints == NULL) {
if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
return 1;
}
/* Create a set of 1D rspl setup points that contains */
/* the residual error */
for (i = 0; i < p->nodp; i++) {
double tv[MXDO]; /* Target output value */
double mv[MXDO]; /* Model output value */
double ev;
xfit_abs_to_rel(p, tv, p->ipoints[i].v);
xfit_shmatsh(p, mv, p->ipoints[i].p);
/* Evaluate the error squared */
if (p->flags & XFIT_FM_INPUT) {
double pp[MXDI];
for (e = 0; e < di; e++)
pp[e] = p->ipoints[i].p[e];
for (f = 0; f < fdi; f++) {
double t1 = tv[f] - mv[f]; /* Error in output */
/* Create input point offset by equivalent delta to output point */
/* error, in proportion to the partial derivatives for that output. */
for (e = 0; e < di; e++)
pp[e] += t1 * p->piv[i].ide[f][e];
}
ev = p->to_de2(p->cntx2, pp, p->ipoints[i].p);
} else {
ev = p->to_de2(p->cntx2, mv, tv);
}
p->rpoints[i].v[0] = ev;
p->rpoints[i].w = p->ipoints[i].w;
}
/* Do each input axis in turn */
for (ee = 0; ee < p->di; ee++) {
rspl *resid;
double imin[1],imax[1],omin[1],omax[1];
int resres[1] = { 1024 };
#define NPGP 100
mcv *posc;
mcvco pgp[NPGP];
double vo, vs;
double *pms;
/* Create a rspl that plots the residual error */
/* vs the axis value */
for (i = 0; i < p->nodp; i++)
p->rpoints[i].p[0] = p->ipoints[i].p[ee];
imin[0] = in_min[ee];
imax[0] = in_max[ee];
omin[0] = 0.0;
omax[0] = 0.0;
if ((resid = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL)
return 1;
resid->fit_rspl_w(resid, RSPLFLAGS, p->rpoints, p->nodp, imin, imax, resres,
omin, omax, 2.0, NULL, NULL);
#ifdef DEBUG_PLOT
{
#define XRES 100
double xx[XRES];
double y1[XRES];
printf("Input residual error channel %d\n",ee);
for (i = 0; i < XRES; i++) {
co pp;
double x;
x = i/(double)(XRES-1);
xx[i] = x = x * (imax[0] - imin[0]) + imin[0];
pp.p[0] = xx[i];
resid->interp(resid, &pp);
y1[i] = pp.v[0];
if (y1[i] < 0.0)
y1[i] = 0.0;
y1[i] = pow(y1[i], 0.5); /* Convert from error^2 to error */
}
do_plot(xx,y1,NULL,NULL,XRES);
}
#endif /* DEBUG_PLOT */
/* Create a set of guide points that contain */
/* the accumulated residual error vs. the axis value */
for (i = 0; i < NPGP; i++) {
co pp;
double vv;
pp.p[0] = i/(NPGP-1.0);
resid->interp(resid, &pp);
pgp[i].p = (pp.p[0] - in_min[ee])/(in_max[ee] - in_min[ee]);
pgp[i].w = 1.0;
vv = pp.v[0];
if (vv < 0.0)
vv = 0.0;
vv = pow(vv, 0.5); /* Convert from error^2 to error */
vv += PSHAPE_MINE; /* In case error is near zero */
vv = pow(vv, PSHAPE_DIST); /* Agressivness of grid distribution */
if (i == 0)
pgp[i].v = vv;
else
pgp[i].v = pgp[i-1].v + vv;
}
resid->del(resid);
/* Normalize the output range */
vo = pgp[0].v;
vs = pgp[NPGP-1].v - vo;
for (i = 0; i < NPGP; i++) {
pgp[i].v = (pgp[i].v - vo)/vs;
/* Apply any dark emphasis */
if (demph > 1.0) {
pgp[i].v = icx_powlike(pgp[i].v, 1.0/demph);
}
}
/* Fit the non-monotonic parameters to the guide points */
if ((posc = new_mcv_noos()) == NULL)
return 1;
posc->fit(posc, 0, p->iluord[ee], pgp, NPGP, 0.1); // ~~99
#ifdef DEBUG_PLOT
{
#define XRES 100
double xx[XRES];
double y1[XRES];
printf("Position curve %d\n",ee);
for (i = 0; i < XRES; i++) {
xx[i] = i/(double)(XRES-1);
y1[i] = posc->interp(posc, xx[i]);
}
do_plot(xx,y1,NULL,NULL,XRES);
}
#endif /* DEBUG_PLOT */
/* Transfer parameters to xfit pos (skip offset and scale) */
posc->get_params(posc, &pms);
for (i = 0; i < p->iluord[ee]; i++) {
p->v[p->pos_offs[ee] + i] = pms[i+2];
}
//p->v[p->in_offs[ee]] = -1.5;
free(pms);
posc->del(posc);
}
}
#ifdef DEBUG
printf("Final parameters:\n");
dump_xfit(p);
#endif
#ifdef DEBUG_PLOT
{
#define XRES 100
double xx[XRES];
double y1[XRES];
for (e = 0; e < p->di; e++) {
printf("Input position curve channel %d\n",e);
for (i = 0; i < XRES; i++) {
double x;
x = i/(double)(XRES-1);
xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
y1[i] = xfit_poscurve(p, x, e);
}
do_plot(xx,y1,NULL,NULL,XRES);
}
for (e = 0; e < p->di; e++) {
printf("Input shape curve channel %d\n",e);
for (i = 0; i < XRES; i++) {
double x;
x = i/(double)(XRES-1);
xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
y1[i] = xfit_shpcurve(p, x, e);
}
do_plot(xx,y1,NULL,NULL,XRES);
}
for (e = 0; e < p->di; e++) {
printf("Combined input curve channel %d\n",e);
for (i = 0; i < XRES; i++) {
double x;
x = i/(double)(XRES-1);
xx[i] = x = x * (p->in_max[e] - p->in_min[e]) + p->in_min[e];
y1[i] = p->incurve(p, x, e);
}
do_plot(xx,y1,NULL,NULL,XRES);
}
for (f = 0; f < p->fdi; f++) {
printf("Output curve channel %d\n",f);
for (i = 0; i < XRES; i++) {
double x;
x = i/(double)(XRES-1);
xx[i] = x = x * (p->out_max[f] - p->out_min[f]) + p->out_min[f];
y1[i] = p->outcurve(p, x, f);
}
do_plot(xx,y1,NULL,NULL,XRES);
}
}
#endif /* DEBUG_PLOT */
/* Create final clut rspl using the established pos/shape/output curves */
/* and white point */
if (flags & XFIT_MAKE_CLUT) {
double *ipos[MXDI];
/* Create an in' -> rout' scattered test point set */
if (p->rpoints == NULL) {
if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
return 1;
}
for (i = 0; i < p->nodp; i++) {
p->rpoints[i].w = p->ipoints[i].w;
xfit_inpscurves(p, p->rpoints[i].p, p->ipoints[i].p);
for (f = 0; f < fdi; f++)
p->rpoints[i].v[f] = p->ipoints[i].v[f];
xfit_abs_to_rel(p, p->rpoints[i].v, p->rpoints[i].v);
xfit_invoutcurves(p, p->rpoints[i].v, p->rpoints[i].v);
if (p->skm) {
/* Look up the skeleton value */
double skval[MXDO];
p->skm->lookup(p->skm, skval, p->ipoints[i].p);
xfit_abs_to_rel(p, skval, skval);
xfit_invoutcurves(p, skval, skval);
//printf("~1 point %d at %f %f %f, targ %f %f %f skm %f %f %f\n", i,p->ipoints[i].p[0],p->ipoints[i].p[1],p->ipoints[i].p[2], p->rpoints[i].v[0],p->rpoints[i].v[1],p->rpoints[i].v[2], skval[0], skval[1], skval[2]);
/* Subtract it from value at this point, */
/* so rspl will fit difference to skeleton model */
for (f = 0; f < fdi; f++)
p->rpoints[i].v[f] -= skval[f];
}
//printf("~1 point %d, w %f, %f %f %f %f -> %f %f %f\n",i,p->rpoints[i].w,p->rpoints[i].p[0], p->rpoints[i].p[1], p->rpoints[i].p[2], p->rpoints[i].p[3],p->rpoints[i].v[0], p->rpoints[i].v[1], p->rpoints[i].v[2]);
}
/* Create ipos[] arrays, that hold the shaper space */
/* grid position due to the positioning curves. */
/* This tells the rspl scattered data interpolator */
/* the grid spacing that smoothness should be */
/* measured against. */
#ifdef DEBUG
printf("~1 about to setup ipos\n");
#endif
for (e = 0; e < p->di; e++) {
//printf("~1 e = %d\n",e);
if ((ipos[e] = (double *)malloc((p->gres[e]) * sizeof(double))) == NULL)
return 1;
//printf("~1 about to do %d spans\n",p->gres[e]);
for (i = 0; i < p->gres[e]; i++) {
double cv;
cv = (double)i/p->gres[e];
//printf("~1 i = %d, pos space = %f\n",i,cv);
/* Inverse lookup grid position through positioning curve */
/* to give device space type value */
cv = icxInvTransFunc(p->v + p->pos_offs[e], p->iluord[e], cv);
//printf("~1 dev space = %f\n",cv);
/* Forward lookup device type value through the shaper curve */
/* to give value in shape linearized space. */
cv = icxTransFunc(p->v + p->shp_offs[e], p->iluord[e], cv);
//printf("~1 shape space = %f\n",cv);
ipos[e][i] = cv;
#ifdef DEBUG
printf("~1 ipos[%d][%d] = %f\n",e,i,cv);
#endif
}
}
if (p->clut != NULL)
p->clut->del(p->clut);
if ((p->clut = new_rspl(RSPL_NOFLAGS, di, fdi)) == NULL)
return 1;
if (p->verb)
printf("Create final clut from scattered data\n");
#ifdef EXTEND_GRID
#define XN EXTEND_GRID_BYN
/* Try increasing the grid by one row all around */
{
#pragma message("!!!!!!!!!!!! Experimental rspl fitting resolution !!!!!!!!!")
double xin_min[MXDI];
double xin_max[MXDI];
int xgres[MXDI];
double del;
double *xipos[MXDI];
for (e = 0; e < p->di; e++) {
del = (in_max[e] - in_min[e])/(gres[e]-1.0); /* Extension */
xin_min[e] = in_min[e] - XN * del;
xin_max[e] = in_max[e] + XN * del;
xgres[e] = gres[e] + 2 * XN;
//printf("~1 xgres %d, gres %d\n",xgres[e], gres[e]);
if ((xipos[e] = (double *)malloc((xgres[e]) * sizeof(double))) == NULL)
return 1;
for (i = 0; i < xgres[e]; i++) {
if (i < XN) { /* Extrapolate bottom */
xipos[e][i] = ipos[e][0] - (XN - i) * (ipos[e][1] - ipos[e][0]);
//printf("~1 xipos[%d] %f from ipos[%d] %f and ipos[%d] %f\n",i,xipos[e][i],0,ipos[e][0],1,ipos[e][1]);
} else if (i >= (xgres[e]-XN)) { /* Extrapolate top */
xipos[e][i] = ipos[e][gres[e]-1] + (i - xgres[e] + XN + 1) * (ipos[e][gres[e]-1] - ipos[e][gres[e]-2]);
//printf("~1 xipos[%d] %f from ipos[%d] %f and ipos[%d] %f\n",i,xipos[e][i],gres[e]-1,ipos[e][gres[e]-1],gres[e]-2,ipos[e][gres[e]-2]);
} else {
xipos[e][i] = ipos[e][i-XN];
//printf("~1 xipos[%d] %f from ipos[%d] %f\n",i,xipos[e][i],i-XN,ipos[e][i-XN]);
}
}
}
p->clut->fit_rspl_w(p->clut, rsplflags, p->rpoints, p->nodp, xin_min, xin_max, xgres,
out_min, out_max, smooth, oavgdev, xipos);
for (e = 0; e < p->di; e++) {
free(xipos[e]);
}
}
#undef XN
#else
// if (p->skm) {
// /* This doesn't seem to work as well as some explicit neutral axis points.. */
// p->clut->fit_rspl_w_df(p->clut, rsplflags, p->rpoints, p->nodp, in_min, in_max, gres,
// out_min, out_max, smooth, oavgdev, ipos, 1.0, (void *)p, skm_weak);
// } else
/* Normal multi-d scattered point fitting */
p->clut->fit_rspl_w(p->clut, rsplflags, p->rpoints, p->nodp, in_min, in_max, gres,
out_min, out_max, smooth, oavgdev, ipos);
#endif
if (p->verb)
printf("\n");
for (e = 0; e < p->di; e++)
free(ipos[e]);
/* If we used a skeleton model, add it into the resulting rspl values */
if (p->skm) {
/* Undo the input point change to allow diagnostic code to work */
for (i = 0; i < p->nodp; i++) {
/* Look up the skeleton value */
double skval[MXDO];
p->skm->lookup(p->skm, skval, p->ipoints[i].p);
xfit_abs_to_rel(p, skval, skval);
xfit_invoutcurves(p, skval, skval);
/* Subtract it from value at this point, */
/* so rspl will fit difference to skeleton model */
for (f = 0; f < fdi; f++)
p->rpoints[i].v[f] += skval[f];
}
/* Undo the skm from the resultant rspl */
p->clut->re_set_rspl(
p->clut, /* this */
0, /* Combination of flags */
(void *)p, /* Opaque function context */
skm_rspl_out /* Function to set from */
);
}
/* The overall device to absolute conversion is now what we want */
/* (as dictated by the points, weighting and best fit), */
/* but we need to adjust the device to relative conversion */
/* to make device white map exactly to D50, without touching */
/* the overall absolute behaviour. */
if (p->flags & XFIT_OUT_WP_REL) {
co wcc; /* device white + aprox rel. white */
icmXYZNumber _wp; /* Uncorrected dw maps to _wp */
if (p->verb)
printf("Doing White point fine tune:\n");
/* See what the relative and absolute white point has turned out to be, */
/* by looking up the device white in the current conversion */
xfit_inpscurves(p, wcc.p, dw);
p->clut->interp(p->clut, &wcc);
xfit_outcurves(p, wcc.v, wcc.v);
if (p->flags & XFIT_OUT_LAB)
icmLab2XYZ(&icmD50, wcc.v, wcc.v);
if (p->verb) {
double labwp[3];
icmXYZ2Lab(&icmD50, labwp, wcc.v);
printf("Before fine tune, rel WP = XYZ %f %f %f, Lab %f %f %f\n",
wcc.v[0], wcc.v[1],wcc.v[2], labwp[0], labwp[1], labwp[2]);
}
/* Matrix needed to correct approx rel wp to target D50 */
icmAry2XYZ(_wp, wcc.v); /* Aprox relative target white point */
if (p->picc != NULL) /* Correction */
p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, NULL, p->cmat, icmD50, _wp);
else
icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->cmat);
/* Compute the actual white point, and return it to caller */
icmMulBy3x3(wp, p->toAbs, wcc.v);
icmAry2Ary(p->wp, wp);
/* Apply correction to fine tune rspl data. */
/* NOTE: this doesn't always give us a perfect D50 white for */
/* Lab PCS input profiles because the dev white may land */
/* within a cell, and the clipping of Lab PCS values in the grid */
/* may introduce errors in the interpolated value. */
p->clut->re_set_rspl(
p->clut, /* this */
0, /* Combination of flags */
(void *)p, /* Opaque function context */
conv_rspl_out /* Function to set from */
);
/* Fix absolute conversions to leave absolute response unchanged. */
icmAry2XYZ(_wp, wp); /* Actual white point */
if (p->picc != NULL) {
p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, p->toAbs, p->fromAbs, icmD50, _wp);
} else {
icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
}
if (p->verb) {
double labwp[3];
/* Lookup white again */
xfit_inpscurves(p, wcc.p, dw);
p->clut->interp(p->clut, &wcc);
xfit_outcurves(p, wcc.v, wcc.v);
if (p->flags & XFIT_OUT_LAB)
icmLab2XYZ(&icmD50, wcc.v, wcc.v);
icmXYZ2Lab(&icmD50, labwp, wcc.v);
printf("After fine tune, rel WP = XYZ %f %f %f, Lab %f %f %f\n",
wcc.v[0], wcc.v[1], wcc.v[2], labwp[0], labwp[1], labwp[2]);
printf(" abs WP = XYZ %s, Lab %s\n", icmPdv(3, wp), icmPLab(wp));
}
}
/* Create default wpscale */
if (wpscale < 0.0) {
wpscale = 1.0;
} else {
if (p->verb) {
printf("White manual point scale %f\n", wpscale);
}
}
/* If we are going to auto scale the WP to avoid clipping */
/* cLUT values above the WP: */
if ((p->flags & XFIT_OUT_WP_REL_US) == XFIT_OUT_WP_REL_US) {
co wcc;
double bw[3];
icmXYZNumber _wp;
double uswpscale = 1.0;
double mxd, mxY;
double ndw[3];
/* See what device space gamut boundary white (ie. 1,1,1) maps to */
xfit_inpscurves(p, wcc.p, dgw);
p->clut->interp(p->clut, &wcc);
xfit_outcurves(p, wcc.v, wcc.v);
if (p->flags & XFIT_OUT_LAB)
icmLab2XYZ(&icmD50, wcc.v, wcc.v);
icmMulBy3x3(wcc.v, p->toAbs, wcc.v); /* Convert to absolute */
mxY = wcc.v[1];
icmCpy3(bw, wcc.v);
//printf("~1 1,1,1 Y = %f\n",wcc.v[1]);
/* See what the device white point value scaled to 1 produces */
mxd = -1.0;
for (e = 0; e < p->di; e++) {
if (dw[e] > mxd)
mxd = dw[e];
}
for (e = 0; e < p->di; e++)
ndw[e] = dw[e]/mxd;
xfit_inpscurves(p, wcc.p, ndw);
p->clut->interp(p->clut, &wcc);
xfit_outcurves(p, wcc.v, wcc.v);
if (p->flags & XFIT_OUT_LAB)
icmLab2XYZ(&icmD50, wcc.v, wcc.v);
icmMulBy3x3(wcc.v, p->toAbs, wcc.v); /* Convert to absolute */
//printf("~1 ndw = %f %f %f Y = %f\n",ndw[0],ndw[1],ndw[2],wcc.v[1]);
if (wcc.v[1] > mxY) {
mxY = wcc.v[1];
icmCpy3(bw, wcc.v);
}
/* Compute WP scale factor needed to fit mxY */
if (mxY > wp[1]) {
uswpscale = mxY/wp[1];
wpscale *= uswpscale;
if (p->verb) {
printf("Dev boundary white XYZ %s, scale WP by %f, total WP scale %f\n",
icmPdv(3, bw), uswpscale, wpscale);
}
}
}
/* If the scaled WP would have Y > 1.0, clip it to 1.0 */
if (p->flags & XFIT_CLIP_WP) {
if ((wp[1] * wpscale) > 1.0) {
wpscale = 1.0/wp[1]; /* Make wp Y = 1.0 */
if (p->verb) {
printf("WP Y would ve > 1.0. scale by %f to clip it\n",wpscale);
}
}
}
/* Apply our total wp scale factor */
if (wpscale != 1.0) {
icmXYZNumber _wp;
/* Create inverse scaling matrix for relative rspl data */
icmSetUnity3x3(p->cmat);
icmScale3x3(p->cmat, p->cmat, 1.0/wpscale);
/* Inverse scale the rspl */
p->clut->re_set_rspl(
p->clut, /* this */
0, /* Combination of flags */
(void *)p, /* Opaque function context */
conv_rspl_out /* Function to set from */
);
/* Scale the WP */
icmScale3(wp, wp, wpscale);
/* return scaled white point to caller */
icmAry2Ary(p->wp, wp);
/* Fix absolute conversions to leave absolute response unchanged. */
icmAry2XYZ(_wp, wp); /* Actual white point */
icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->fromAbs);
icmChromAdaptMatrix(ICM_CAM_BRADFORD, _wp, icmD50, p->toAbs);
}
/* Clip any values in the grid over D50 L to D50 */
if ((p->flags & XFIT_OUT_WP_REL_C) == XFIT_OUT_WP_REL_C) {
/* Compute the rspl D50 value to avoid calc in clip_rspl_out() */
if (p->flags & XFIT_OUT_LAB) {
p->cmat[0][0] = 100.0;
p->cmat[0][1] = 0.0;
p->cmat[0][2] = 0.0;
} else {
icmXYZ2Ary(p->cmat[0], icmD50);
}
xfit_invoutcurves(p, p->cmat[0], p->cmat[0]);
if (p->verb)
printf("Clipping any cLUT grid points with Y > 1 to D50\n");
p->clut->re_set_rspl(
p->clut, /* this */
0, /* Combination of flags */
(void *)p, /* Opaque function context */
clip_rspl_out /* Function to set from */
);
}
/* Force black point to given value */
// if (bpo != NULL) {
// co tv;
// int rv;
//
// xfit_inpscurves(p, tv.p, bpo->p);
//
// xfit_XYZ_abs_to_rel(p, tv.v, bpo->v);
// xfit_invoutcurves(p, tv.v, tv.v);
//printf("~1 xfit: fine after curves black at %f %f %f to %f %f %f\n",
//tv.p[0], tv.p[1], tv.p[2], tv.v[0], tv.v[1], tv.v[2]);
// rv = p->clut->tune_value(p->clut, &tv);
// if (rv != 0)
// warning("Black Point Override failed - clipping");
// }
#ifdef SPECIAL_FORCE
/* Replace the rspl nodes with ones directly computed */
/* from the synthetic linear RGB->XYZ model */
{
double twp[3];
icmXYZNumber _wp;
/* See what current device white maps to */
twp[0] = twp[1] = twp[2] = 1.0;
domodel(twp, twp);
icmAry2XYZ(_wp, twp);
/* Matrix needed to correct to D50 */
icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, _wp, p->cmat);
p->clut->re_set_rspl(
p->clut, /* this */
0, /* Combination of flags */
(void *)p, /* Opaque function context */
set_rspl_out1 /* Function to set from */
);
}
#endif /* SPECIAL_FORCE */
/* Evaluate the residual error now, with the rspl in place */
#ifdef DEBUG_PLOT
{
int ee;
double maxe = 0.0;
double avee = 0.0;
/* Allocate in->rout duplicate point set */
if (p->rpoints == NULL) {
if ((p->rpoints = (cow *)malloc(p->nodp * sizeof(cow))) == NULL)
return 1;
}
/* Create a set of 1D rspl setup points that contains */
/* the residual error */
for (i = 0; i < p->nodp; i++) {
double tv[MXDO]; /* Target output value */
double mv[MXDO]; /* Model output value */
co pp;
double ev;
xfit_abs_to_rel(p, tv, p->ipoints[i].v);
for (e = 0; e < p->di; e++)
pp.p[e] = p->ipoints[i].p[e];
xfit_inpscurves(p, pp.p, pp.p);
p->clut->interp(p->clut, &pp);
xfit_outcurves(p, pp.v, pp.v);
for (f = 0; f < p->fdi; f++)
mv[f] = pp.v[f];
/* Evaluate the residual error suqared */
if (p->flags & XFIT_FM_INPUT) {
double pp[MXDI];
for (e = 0; e < di; e++)
pp[e] = p->ipoints[i].p[e];
for (f = 0; f < fdi; f++) {
double t1 = tv[f] - mv[f]; /* Error in output */
/* Create input point offset by equivalent delta to output point */
/* error, in proportion to the partial derivatives for that output. */
for (e = 0; e < di; e++)
pp[e] += t1 * p->piv[i].ide[f][e];
}
ev = p->to_de2(p->cntx2, pp, p->ipoints[i].p);
} else {
ev = p->to_de2(p->cntx2, mv, tv);
}
//printf("~1 point %d, loc %f %f %f %f\n",i, p->rpoints[i].p[0], p->rpoints[i].p[1], p->rpoints[i].p[2], p->rpoints[i].p[3]);
//printf(" targ %f %f %f, is %f %f %f\n", tv[0], tv[1], tv[2], mv[0], mv[1], mv[2]);
p->rpoints[i].v[0] = ev;
p->rpoints[i].w = p->ipoints[i].w;
ev = sqrt(ev);
if (ev > maxe)
maxe = ev;
avee += ev;
}
printf("Max resid err = %f, avg err = %f\n",maxe, avee/(double)p->nodp);
/* Evaluate each input axis in turn */
for (ee = 0; ee < p->di; ee++) {
rspl *resid;
double imin[1],imax[1],omin[1],omax[1];
int resres[1] = { 1024 };
/* Create a rspl that gives the residual error squared */
/* vs. the axis value */
for (i = 0; i < p->nodp; i++)
p->rpoints[i].p[0] = p->ipoints[i].p[ee];
imin[0] = in_min[ee];
imax[0] = in_max[ee];
omin[0] = 0.0;
omax[0] = 0.0;
if ((resid = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL)
return 1;
resid->fit_rspl_w(resid, RSPLFLAGS, p->rpoints, p->nodp, imin, imax, resres,
omin, omax, 2.0, NULL, NULL);
{
#define XRES 100
double xx[XRES];
double y1[XRES];
printf("Finale residual error vs. input channel %d\n",ee);
for (i = 0; i < XRES; i++) {
co pp;
double x;
x = i/(double)(XRES-1);
xx[i] = x = x * (imax[0] - imin[0]) + imin[0];
pp.p[0] = xx[i];
resid->interp(resid, &pp);
y1[i] = sqrt(fabs(pp.v[0]));
}
do_plot(xx,y1,NULL,NULL,XRES);
}
resid->del(resid);
}
}
#endif /* DEBUG_PLOT */
/* Special test code to figure out what's wrong with position curves */
#ifdef NEVER
{
double rgb[3], xyz[3];
co pp;
extern int rspldb;
db = 1;
rspldb = 1;
printf("~1 gres = %d %d %d\n", gres[0], gres[1], gres[2]);
for (i = 0; i < 3; i++) {
if (i == 0) {
/* Test point on a grid point */
printf("\n~1 ##########################################\n");
printf("~1 testing input at first diagonal grid point\n");
rgb[0] = 1.0/(gres[0]-1.0);
rgb[1] = 1.0/(gres[0]-1.0);
rgb[2] = 1.0/(gres[0]-1.0);
printf("~1 target rgb' = %f %f %f\n", rgb[0], rgb[1], rgb[2]);
xfit_invinpscurves(p, rgb, rgb);
} else if (i == 1) {
/* Test point half way through a grid point */
printf("\n~1 ##########################################\n");
printf("\n~1 testing half way through diagonal grid point\n");
rgb[0] = 0.5/(gres[0]-1.0);
rgb[1] = 0.5/(gres[0]-1.0);
rgb[2] = 0.5/(gres[0]-1.0);
printf("~1 target rgb' = %f %f %f\n", rgb[0], rgb[1], rgb[2]);
xfit_invinpscurves(p, rgb, rgb);
} else {
printf("\n~1 ##########################################\n");
printf("\n~1 testing worst case point\n");
rgb[0] = 0.039915;
rgb[1] = 0.053148;
rgb[2] = 0.230610;
}
pp.p[0] = rgb[0];
pp.p[1] = rgb[1];
pp.p[2] = rgb[2];
printf("~1 rgb = %f %f %f\n", pp.p[0], pp.p[1], pp.p[2]);
xfit_inpscurves(p, pp.p, pp.p);
printf("~1 rgb' = %f %f %f\n", pp.p[0], pp.p[1], pp.p[2]);
p->clut->interp(p->clut, &pp);
printf("~1 xyz' = %f %f %f\n", pp.v[0], pp.v[1], pp.v[2]);
xfit_outcurves(p, pp.v, pp.v);
printf("~1 xyz = %f %f %f\n", pp.v[0], pp.v[1], pp.v[2]);
/* Synthetic linear rgb->XYZ model */
domodel(xyz, rgb);
/* Apply abs->rel white point adjustment */
icmMulBy3x3(xyz, p->mat, xyz);
printf("~1 ref = %f %f %f, de = %f\n", xyz[0], xyz[1], xyz[2],icmXYZLabDE(&icmD50,xyz,pp.v));
}
db = 0;
rspldb = 0;
// exit(0);
}
#endif /* NEVER */
}
return 0;
}
/* We're done with an xfit */
static void xfit_del(xfit *p) {
if (p->v != NULL)
free(p->v);
if (p->wv != NULL)
free(p->wv);
if (p->sa != NULL)
free(p->sa);
if (p->rpoints != NULL)
free(p->rpoints);
if (p->piv != NULL)
free(p->piv);
if (p->uerrv != NULL)
free(p->uerrv);
free(p);
}
/* Create a transform fitting object */
/* return NULL on error */
xfit *new_xfit(
icc *picc /* ICC profile used to set cone space matrix, NULL for Bradford. */
) {
xfit *p;
if ((p = (xfit *)calloc(1, sizeof(xfit))) == NULL) {
return NULL;
}
p->picc = picc;
/* Set method pointers */
p->fit = xfit_fit;
p->incurve = xfit_inpscurve;
p->invincurve = xfit_invinpscurve;
p->outcurve = xfit_outcurve;
p->invoutcurve = xfit_invoutcurve;
p->del = xfit_del;
return p;
}
Reference in New Issue View Git Blame Copy Permalink