3151 lines
90 KiB
C
3151 lines
90 KiB
C
|
|
|
|
/*
|
|
* 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;
|
|
}
|
|
|
|
|
|
|