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

2104 lines
59 KiB
C
Raw Permalink Blame History

/*
* International Color Consortium color transform expanded support
*
* Author: Graeme W. Gill
* Date: 2/7/00
* Version: 1.00
*
* Copyright 2000, 2003 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 old iccXfm class.
*/
/*
* This module provides the expands icclib functionality
* for matrix profiles.
* This file is #included in xicc.c, to keep its functions private.
*/
/*
* TTBD:
*
* Some of the error handling is crude. Shouldn't use
* error(), should return status.
*
* Should allow for offset in curves - this will greatly improve
* profile quality on non-calibrated displays. See spectro/dispcal.c
* spectro/moncurve.c. Use conjgrad() instead of powell() to speed things up.
* Note that if curves have scale, the scale will have to be
* normalized back to zero by scaling the matrix before storing
* the result in the ICC profile.
*
* Would be nice to be able to specify a specific patch
* as the white one rather than using heuristic to identify it,
* since some pathalogical cases don't work.
*/
#define USE_CIE94_DE /* Use CIE94 delta E measure when creating fit */
/* Weights in shaper parameters, to minimise unconstrained "wiggles" */
#define MXNORDERS 30 /* Maximum shaper harmonic orders to use */
#define XSHAPE_MAG 1.0 /* Overall shaper parameter magnitide */
#define XSHAPE_OFFG 0.1 /* Input offset weights when ord 0 is gamma */
#define XSHAPE_OFFS 1.0 /* Input offset weights when ord 0 is shaper */
#define XSHAPE_HW01 0.01 /* 0 & 1 harmonic weights */
#define XSHAPE_HBREAK 3 /* Harmonic that has HWBR */
#define XSHAPE_HWBR 0.5 /* Base weight of harmonics HBREAK up */
#define XSHAPE_HWINC 0.5 /* Increase in weight for each harmonic above HWBR */
#define XSHAPE_GAMTHR 0.01 /* Input threshold for linear slope below gamma power */
#undef DEBUG /* [und] Extra printfs */
#undef DEBUG_PLOT /* [und] Plot curves */
#undef DEBUG_SPEC /* [und] Debug some specific cases */
#define G_DEB 0 /* [und] g_deb default value */
/* Debug some specific cases (fwd_relpcs_outpcs, bwd_outpcs_relpcs) */
#ifdef DEBUG_SPEC
# undef DBS
# ifdef CNDTRACE
# define DBS(xxx) if (p->trace) printf xxx ;
# else
# define DBS(xxx) printf xxx ;
# endif
#else
# undef DBS
# define DBS(xxx)
#endif
/* ========================================================= */
/* Forward and Backward Matrix type conversion */
/* Return 0 on success, 1 if clipping occured, 2 on other error */
/* Individual components of Fwd conversion: */
/* Because icm_lu4 sets up component conversions in the requested direction, */
/* but icxLuMatrix expects component conversions as if the the matrix is */
/* fwd, we have if (dir) code to reverse things here. */
static int
icxLuMatrixFwd_curve (
icxLuMatrix *p, /* This */
double *out, /* Vector of output values */
double *in /* Vector of input values */
) {
icmPe_lurv prv;
if (p->dir) {
prv = p->plu->output_fmt_bwd(p->plu, out, in); /* Should be NOP */
prv |= p->plu->output_pch_bwd(p->plu, out, out);
} else {
prv = p->plu->input_fmt_fwd(p->plu, out, in); /* Should be NOP */
prv |= p->plu->input_pch_fwd(p->plu, out, out);
}
return LUE2XLUE(prv);
}
static int
icxLuMatrixFwd_matrix (
icxLuMatrix *p, /* This */
double *out, /* Vector of output values */
double *in /* Vector of input values */
) {
icmPe_lurv prv;
if (p->dir)
prv = p->plu->core5_bwd(p->plu, out, in);
else
prv = p->plu->core5_fwd(p->plu, out, in);
return LUE2XLUE(prv);
}
static int
icxLuMatrixFwd_abs (
icxLuMatrix *p, /* This */
double *out, /* Vector of output values */
double *in /* Vector of input values */
) {
int rv = 0;
icmPe_lurv prv;
if (p->dir) {
prv = p->plu->input_pch_bwd(p->plu, out, in); /* Should be NOP */
prv |= p->plu->input_fmt_bwd(p->plu, out, out);
} else {
prv = p->plu->output_pch_fwd(p->plu, out, in); /* Should be NOP */
prv |= p->plu->output_fmt_fwd(p->plu, out, out);
}
rv = LUE2XLUE(prv);
if (p->pcs == icxSigJabData)
p->cam->XYZ_to_cam(p->cam, out, out);
return rv;
}
/* Overall Fwd conversion */
/* Note that the overall conversion is in the requested direction, */
/* as the setup code swaps icxLuMatrixFwd_lookup/icxLuMatrixBwd_lookup */
/* as needed. */
static int
icxLuMatrixFwd_lookup (
icxLuBase *pp, /* This */
double *out, /* Vector of output values */
double *in /* Vector of input values */
) {
int rv = 0;
icxLuMatrix *p = (icxLuMatrix *)pp;
rv |= icxLuMatrixFwd_curve(p, out, in);
rv |= icxLuMatrixFwd_matrix(p, out, out);
rv |= icxLuMatrixFwd_abs(p, out, out);
return rv;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Given a relative XYZ or Lab PCS value, convert in the fwd direction into */
/* the nominated output PCS (ie. Absolute, Jab etc.) */
/* (This is used in generating gamut compression in B2A tables) */
void icxLuMatrix_fwd_relpcs_outpcs(
icxLuBase *pp,
icColorSpaceSignature is, /* Input space, XYZ or Lab */
double *out, double *in) {
icxLuMatrix *p = (icxLuMatrix *)pp;
if (is == icSigLabData && p->natpcs == icSigXYZData) {
DBS(("Matrix_fwd_relpcs_outpcs: Lab in = %s\n", icmPdv(p->inputChan, in)));
icmLab2XYZ(&icmD50, out, in);
DBS(("Matrix_fwd_relpcs_outpcs: XYZ = %s\n", icmPdv(p->outputChan, out)));
icxLuMatrixFwd_abs(p, out, out);
DBS(("Matrix_fwd_relpcs_outpcs: abs XYZ = %s\n", icmPdv(p->outputChan, out)));
} else if (is == icSigXYZData && p->natpcs == icSigLabData) {
DBS(("Matrix_fwd_relpcs_outpcs: XYZ in = %s\n", icmPdv(p->inputChan, in)));
icmXYZ2Lab(&icmD50, out, in);
DBS(("Matrix_fwd_relpcs_outpcs: Lab = %s\n", icmPdv(p->outputChan, out)));
icxLuMatrixFwd_abs(p, out, out);
DBS(("Matrix_fwd_relpcs_outpcs: abs Lab = %s\n", icmPdv(p->outputChan, out)));
} else {
DBS(("Matrix_fwd_relpcs_outpcs: rel = %s\n", icmPdv(p->inputChan, in)));
icxLuMatrixFwd_abs(p, out, in);
DBS(("Matrix_fwd_relpcs_outpcs: abs = %s\n", icmPdv(p->outputChan, out)));
}
}
/* - - - - - - - - - - - - - - - - - - - - - */
/* Individual components of Bwd conversion: */
static int
icxLuMatrixBwd_abs (
icxLuMatrix *p, /* This */
double *out, /* Vector of output values */
double *in /* Vector of input values */
) {
icmPe_lurv prv;
int rv = 0;
if (p->pcs == icxSigJabData) {
p->cam->cam_to_XYZ(p->cam, out, in);
/* Hack to prevent CAM02 weirdness being amplified by */
/* any later per channel clipping. */
/* Limit -Y to non-stupid values by scaling */
if (out[1] < -0.1) {
out[0] *= -0.1/out[1];
out[2] *= -0.1/out[1];
out[1] = -0.1;
}
if (p->dir) {
prv = p->plu->input_fmt_fwd(p->plu, out, out);
prv |= p->plu->input_pch_fwd(p->plu, out, out); /* Should be NOP */
} else {
prv = p->plu->output_fmt_bwd(p->plu, out, out);
prv |= p->plu->output_pch_bwd(p->plu, out, out); /* Should be NOP */
}
rv = LUE2XLUE(prv);
} else {
if (p->dir) {
prv = p->plu->input_fmt_fwd(p->plu, out, in);
prv |= p->plu->input_pch_fwd(p->plu, out, out); /* Should be NOP */
} else {
prv = p->plu->output_fmt_bwd(p->plu, out, in);
prv |= p->plu->output_pch_bwd(p->plu, out, out); /* Should be NOP */
}
rv = LUE2XLUE(prv);
}
return rv;
}
static int
icxLuMatrixBwd_matrix (
icxLuMatrix *p, /* This */
double *out, /* Vector of output values */
double *in /* Vector of input values */
) {
icmPe_lurv prv;
if (p->dir)
prv = p->plu->core5_fwd(p->plu, out, in);
else
prv = p->plu->core5_bwd(p->plu, out, in);
return LUE2XLUE(prv);
}
static int
icxLuMatrixBwd_curve (
icxLuMatrix *p, /* This */
double *out, /* Vector of output values */
double *in /* Vector of input values */
) {
icmPe_lurv prv;
if (p->dir) {
prv = p->plu->output_pch_fwd(p->plu, out, in);
prv |= p->plu->output_fmt_fwd(p->plu, out, out); /* Should be NOP */
} else {
prv = p->plu->input_pch_bwd(p->plu, out, in);
prv |= p->plu->input_fmt_bwd(p->plu, out, out); /* Should be NOP */
}
return LUE2XLUE(prv);
}
/* Overall Bwd conversion */
/* Note that the overall conversion is in the requested direction, */
/* as the setup code swaps icxLuMatrixFwd_lookup/icxLuMatrixBwd_lookup */
/* as needed. */
static int
icxLuMatrixBwd_lookup (
icxLuBase *pp, /* This */
double *out, /* Vector of output values */
double *in /* Vector of input values */
) {
int rv = 0;
icxLuMatrix *p = (icxLuMatrix *)pp;
rv |= icxLuMatrixBwd_abs(p, out, in);
rv |= icxLuMatrixBwd_matrix(p, out, out);
rv |= icxLuMatrixBwd_curve(p, out, out);
return rv;
}
static void
icxLuMatrix_free(
icxLuBase *p
) {
p->plu->del(p->plu);
if (p->cam != NULL)
p->cam->del(p->cam);
free(p);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Given a nominated output PCS (ie. Absolute, Jab etc.), convert it in the bwd */
/* direction into a relative XYZ or Lab PCS value */
/* (This is used in generating gamut compression in B2A tables) */
void icxLuMatrix_bwd_outpcs_relpcs(
icxLuBase *pp,
icColorSpaceSignature os, /* Output space, XYZ or Lab */
double *out, double *in) {
icxLuMatrix *p = (icxLuMatrix *)pp;
DBS(("Matrix_bwd_outpcs_relpcs: rel = %s\n", icmPdv(p->outputChan, in)));
icxLuMatrixBwd_abs(p, out, in);
DBS(("Matrix_bwd_outpcs_relpcs: abs = %s\n", icmPdv(p->inputChan, out)));
if (os == icSigXYZData && p->natpcs == icSigLabData) {
icmLab2XYZ(&icmD50, out, out);
DBS(("Matrix_bwd_outpcs_relpcs: XYZ = %s\n", icmPdv(p->inputChan, out)));
} else if (os == icSigLabData && p->natpcs == icSigXYZData) {
icmXYZ2Lab(&icmD50, out, out);
DBS(("Matrix_bwd_outpcs_relpcs: Lab = %s\n", icmPdv(p->inputChan, out)));
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
static gamut *icxLuMatrixGamut(icxLuBase *plu, double detail);
static icxCuspMap *icxLuMatrixCuspMap(icxLuBase *plu, int res) { return NULL; };
/* Do the basic icxLuMatrix creation and initialisation */
static icxLuMatrix *
alloc_icxLuMatrix(
xicc *xicp,
icmLuSpace *plu, /* Pointer to Lu we are expanding (ours) */
int dir, /* 0 = fwd, 1 = bwd */
int flags /* clip, merge flags */
) {
icxLuMatrix *p;
if ((p = (icxLuMatrix *) calloc(1,sizeof(icxLuMatrix))) == NULL)
return NULL;
p->lutype = icxLuMatrixType;
p->pp = xicp;
p->plu = plu;
p->del = icxLuMatrix_free;
p->lutspaces = icxLutSpaces;
p->spaces = icxLuSpaces;
p->get_native_ranges = icxLu_get_native_ranges;
p->get_ranges = icxLu_get_ranges;
p->efv_wh_bk_points = icxLuEfv_wh_bk_points;
p->get_gamut = icxLuMatrixGamut;
p->get_cuspmap = icxLuMatrixCuspMap;
p->fwd_relpcs_outpcs = icxLuMatrix_fwd_relpcs_outpcs;
p->bwd_outpcs_relpcs = icxLuMatrix_bwd_outpcs_relpcs;
p->nearclip = 0; /* Set flag defaults */
p->mergeclut = 0;
p->intsep = 0;
p->dir = dir;
p->fwd_lookup = icxLuMatrixFwd_lookup;
p->fwd_curve = icxLuMatrixFwd_curve;
p->fwd_matrix = icxLuMatrixFwd_matrix;
p->fwd_abs = icxLuMatrixFwd_abs;
p->bwd_lookup = icxLuMatrixBwd_lookup;
p->bwd_abs = icxLuMatrixBwd_abs;
p->bwd_matrix = icxLuMatrixBwd_matrix;
p->bwd_curve = icxLuMatrixBwd_curve;
if (dir) { /* Bwd */
p->lookup = icxLuMatrixBwd_lookup;
p->inv_lookup = icxLuMatrixFwd_lookup;
} else {
p->lookup = icxLuMatrixFwd_lookup;
p->inv_lookup = icxLuMatrixBwd_lookup;
}
/* There are no matrix specific flags */
p->flags = flags;
return p;
}
/* We setup fwd and bwd component conversions as if this was a fwd conversion, */
/* but setup the asked for overal conversion in the resquested direction. */
static icxLuBase *
new_icxLuMatrix(
xicc *xicp,
int flags, /* clip, merge flags */
icmLuSpace *plu, /* Pointer to Lu we are expanding */
icmLookupFunc func, /* Functionality requested */
icRenderingIntent intent, /* Rendering intent */
icColorSpaceSignature pcsor, /* PCS override (0 = def) */
icxViewCond *vc, /* Viewing Condition (NULL if pcsor is not CIECAM) */
int dir /* 0 = fwd, 1 = bwd */
) {
icxLuMatrix *p;
// ~8 should merge alloc_icxLuMatrix() into new_icxLuMatrix()
/* Do basic creation and initialisation */
if ((p = alloc_icxLuMatrix(xicp, plu, dir, flags)) == NULL)
return NULL;
p->func = func;
/* Init the CAM model */
if (pcsor == icxSigJabData) {
if (vc != NULL) /* One has been provided */
p->vc = *vc; /* Copy the structure */
else
xicc_enum_viewcond(xicp, &p->vc, -1, NULL, 0, NULL); /* Use a default */
p->cam = new_icxcam(cam_default);
p->cam->set_view_vc(p->cam, &p->vc);
} else {
p->cam = NULL;
}
/* Remember the effective intent */
p->intent = intent;
/* In general the native and effective space info of the icx will be the same as the */
/* underlying icm lookup object. */
{
icmCSInfo ini, outi, pcsi;
/* Get details of internal, native color space */
plu->native_spaces(p->plu, &ini, &outi, &pcsi);
p->natis = ini.sig;
p->natos = outi.sig;
p->natpcs = pcsi.sig;
icmCpyN(p->ninmin, ini.min, ini.nch);
icmCpyN(p->ninmax, ini.max, ini.nch);
icmCpyN(p->noutmin, outi.min, outi.nch);
icmCpyN(p->noutmax, outi.max, outi.nch);
/* Get details of external, effective color spaces */
plu->spaces(p->plu, &ini, &outi, &pcsi, NULL, NULL, NULL, NULL, NULL, NULL);
p->ins = ini.sig;
p->inputChan = ini.nch;
p->outs = outi.sig;
p->outputChan = outi.nch;
p->pcs = pcsi.sig;
icmCpyN(p->inmin, ini.min, ini.nch);
icmCpyN(p->inmax, ini.max, ini.nch);
icmCpyN(p->outmin, outi.min, outi.nch);
icmCpyN(p->outmax, outi.max, outi.nch);
}
/* Override with pcsor */
if (pcsor == icxSigJabData) {
p->pcs = pcsor;
if (func == icmBwd || func == icmGamut || func == icmPreview)
p->ins = pcsor;
if (func == icmFwd || func == icmPreview)
p->outs = pcsor;
}
/* If we have a Jab PCS override, reflect this in the effective icx range. */
/* Note that the ab ranges are nominal. They will exceed this range */
/* for colors representable in L*a*b* PCS */
if (p->ins == icxSigJabData) {
p->inmin[0] = 0.0; p->inmax[0] = 100.0;
p->inmin[1] = -128.0; p->inmax[1] = 128.0;
p->inmin[2] = -128.0; p->inmax[2] = 128.0;
} else if (p->outs == icxSigJabData) {
p->outmin[0] = 0.0; p->outmax[0] = 100.0;
p->outmin[1] = -128.0; p->outmax[1] = 128.0;
p->outmin[2] = -128.0; p->outmax[2] = 128.0;
}
return (icxLuBase *)p;
}
/* ========================================================== */
/* xicc creation code */
/* ========================================================== */
/* Context for figuring input curves */
typedef struct {
rspl *r; /* Device -> PCS rspl */
int linear; /* Flag */
double nmin, nmax; /* PCS End points to linearise */
double min, max; /* device End points to linearise */
} mxinctx;
#define NPARMS (9 + 6 + 3 * MXNORDERS)
/* Context for optimising matrix */
typedef struct {
int verb; /* Verbose */
int optdim; /* Optimisation dimensions */
int isLinear; /* NZ if no curves, fixed Gamma = 1.0 */
int isGamma; /* NZ if gamma + matrix, else shaper */
int isShTRC; /* NZ if shared TRC */
int shape0gam; /* NZ if zero'th order shaper should be gamma function */
int norders; /* Number of shaper orders */
int clipbw; /* Prevent white > 1 and -ve black */
int clipprims; /* Prevent primaries going -ve */
double smooth; /* Shaper smoothing factor (nominal = 1.0) */
double dscale; /* Scale device values */
double v[NPARMS]; /* Holder for parameters */
double sa[NPARMS]; /* Initial search area */
/* Rest are matrix : */
/* 0 1 2 R X */
/* 3 4 5 * G = Y */
/* 6 7 8 B Z */
/* For Gamma: */
/* 9, 10, 11 are gamma */
/* Else for shaper only: */
/* 9, 10, 11 are Input Offset */
/* 12, 13, 14 are Output Offset */
/* 15, 16, 17 are Gamma or 0th harmonics */
/* 18, 19, 20 are 1st harmonics */
/* 21, 22, 23 are 2nd harmonics */
/* 24, 25, 26 etc. */
/* For isShTRC there is only one gamma or set of offsets & harmonics */
icmXYZNumber wp; /* Assumed white point for Lab conversion */
cow *points; /* List of test points as dev->Lab */
int nodp; /* Number of data points */
} mxopt;
/* Per chanel function being optimised */
static void mxmfunc1(mxopt *p, int j, double *v, double *out, double *in) {
double vv, g;
int ps = 3; /* Parameter spacing */
vv = *in * p->dscale;
if (p->isShTRC) {
j = 0;
ps = 1; /* Only one channel */
}
if (p->isLinear) { /* No per channel curve */
*out = vv;
return;
}
if (p->isGamma) { /* Pure Gamma */
/* Apply gamma */
g = v[9 + j];
if (g <= 0.0)
vv = 1.0;
else {
if (vv >= 0.0) {
vv = pow(vv, g);
} else {
vv = -pow(-vv, g);
}
}
} else { /* Add extra shaper parameters */
int ord;
if (p->shape0gam) {
/* Apply input offset */
g = v[9 + j]; /* Offset value */
if (g >= 1.0) {
vv = 1.0;
} else {
vv = g + ((1.0 - g) * vv);
}
/* Apply gamma as order 0 */
g = v[9 + 2 * ps + j];
if (g <= 0.0)
vv = 1.0;
else {
/* Power with straight line at small values */
if (vv >= XSHAPE_GAMTHR) {
vv = pow(vv, g);
} else {
double slope, oth;
oth = pow(XSHAPE_GAMTHR, g); /* Output at input threshold */
slope = g * pow(XSHAPE_GAMTHR, g - 1.0); /* Slope at input threshold */
vv = oth + (vv - XSHAPE_GAMTHR) * slope; /* Straight line */
}
}
}
/* Process all the shaper orders from high to low. */
/* [These shapers were inspired by a Graphics Gem idea */
/* (Gems IV, VI.3, "Fast Alternatives to Perlin's Bias and */
/* Gain Functions, pp 401). */
/* They have the nice properties that they are smooth, and */
/* can't be non-monotonic. The control parameter has been */
/* altered to have a range from -oo to +oo rather than 0.0 to 1.0 */
/* so that the search space is less non-linear. ] */
if (p->shape0gam)
ord = 1;
else
ord = 0;
for (; ord < p->norders; ord++)
{
int nsec; /* Number of sections */
double sec; /* Section */
g = v[9 + 2 * ps + ord * ps + j]; /* Parameter */
nsec = ord + 1; /* Increase sections for each order */
vv *= (double)nsec;
sec = floor(vv);
if (((int)sec) & 1)
g = -g; /* Alternate action in each section */
vv -= sec;
if (g >= 0.0) {
vv = vv/(g - g * vv + 1.0);
} else {
vv = (vv - g * vv)/(1.0 - g * vv);
}
vv += sec;
vv /= (double)nsec;
}
/* (For extrapolation it helps to pin 0 & 1) */
if (p->shape0gam) {
/* Apply output offset */
g = v[9 + 1 * ps + j]; /* Offset value */
if (g >= 1.0) {
vv = 1.0;
} else if (g > 0.0) {
vv = g + ((1.0 - g) * vv);
}
}
}
*out = vv;
}
/* Function being optimised */
static void mxmfunc(mxopt *p, double *v, double *xyz, double *in) {
int j;
double rgb[3];
/* Apply per channel processing */
for (j = 0; j < 3; j++)
mxmfunc1(p, j, v, &rgb[j], &in[j]);
/* Apply matrix */
xyz[0] = v[0] * rgb[0] + v[1] * rgb[1] + v[2] * rgb[2];
xyz[1] = v[3] * rgb[0] + v[4] * rgb[1] + v[5] * rgb[2];
xyz[2] = v[6] * rgb[0] + v[7] * rgb[1] + v[8] * rgb[2];
}
/* return the sum of the squares of the current shaper parameters */
static double xshapmag(
mxopt *p, /* Base of optimisation structure */
double *v /* Pointer to parameters */
) {
double tt, w, tparam = 0.0;
int f, g;
if (p->isGamma) { /* Pure Gamma only */
return 0.0;
}
if (p->isShTRC) {
/* Input offset value */
if (p->shape0gam)
w = XSHAPE_OFFG;
else
w = XSHAPE_OFFS;
tt = v[9];
tt *= tt;
tparam += w * tt;
/* Output offset value */
tt = v[10];
tt *= tt;
tparam += w * tt;
/* Shaper values */
for (f = 0; f < p->norders; f++) {
tt = v[11 + f];
if (f == 0 && p->shape0gam)
tt -= 1.0; /* default is linear */
/* Weigh to suppress ripples */
if (f <= 1) { /* Use XSHAPE_HW01 */
w = XSHAPE_HW01;
} else if (f <= XSHAPE_HBREAK) { /* Blend from XSHAPE_HW01 to XSHAPE_HWBR * smooth */
double bl = (f - 1.0)/(XSHAPE_HBREAK - 1.0);
w = (1.0 - bl) * XSHAPE_HW01 + bl * XSHAPE_HWBR * p->smooth;
} else { /* Use XSHAPE_HWBR * smooth */
w = XSHAPE_HWBR + (f-XSHAPE_HBREAK) * XSHAPE_HWINC * p->smooth;
}
tt *= tt;
tparam += w * tt;
}
return XSHAPE_MAG * tparam;
}
/* Input offset value */
if (p->shape0gam)
w = XSHAPE_OFFG;
else
w = XSHAPE_OFFS;
for (g = 0; g < 3; g++) {
tt = v[9 + g];
tt *= tt;
tparam += w * tt;
}
/* Output ffset value */
for (g = 0; g < 3; g++) {
tt = v[12 + g];
tt *= tt;
tparam += w * tt;
}
/* Shaper values */
for (f = 0; f < p->norders; f++) {
/* Weigh to suppress ripples */
if (f <= 1) {
w = XSHAPE_HW01;
} else if (f <= XSHAPE_HBREAK) {
double bl = (f - 1.0)/(XSHAPE_HBREAK - 1.0);
w = (1.0 - bl) * XSHAPE_HW01 + bl * XSHAPE_HWBR * p->smooth;
} else {
w = XSHAPE_HWBR + (f-XSHAPE_HBREAK) * XSHAPE_HWINC * p->smooth;
}
for (g = 0; g < 3; g++) {
tt = v[15 + 3 * f + g];
if (f == 0 && p->shape0gam)
tt -= 1.0; /* default is linear */
tt *= tt;
tparam += w * tt;
}
}
return XSHAPE_MAG * tparam/3.0;
}
int g_deb = G_DEB;
/* Matrix optimisation function handed to powell() */
static double mxoptfunc(void *edata, double *v) {
mxopt *p = (mxopt *)edata;
double err = 0.0, rv = 0.0, smv;
double xyz[3], lab[3];
int i;
if (g_deb) printf("\n");
for (i = 0; i < p->nodp; i++) {
/* Apply our function */
//printf("%f %f %f -> %f %f %f\n", p->points[i].p[0], p->points[i].p[1], p->points[i].p[2], xyz[0], xyz[1], xyz[2]);
mxmfunc(p, v, xyz, p->points[i].p);
/* Convert to Lab */
icmXYZ2Lab(&p->wp, lab, xyz);
if (g_deb) printf("%d: %f %f %f -> %f %f %f, target %f %f %f, w %f\n", i, p->points[i].p[0], p->points[i].p[1], p->points[i].p[2], lab[0], lab[1], lab[2], p->points[i].v[0], p->points[i].v[1], p->points[i].v[2],p->points[i].w);
/* Accumulate total delta E squared */
#ifdef USE_CIE94_DE
rv += p->points[i].w * icmCIE94sq(lab, p->points[i].v);
#else
rv += p->points[i].w * icmLabDEsq(lab, p->points[i].v);
#endif
}
/* Normalise error to be an average delta E squared */
rv /= (double)p->nodp;
/* Sum with shaper parameters squared, to */
/* minimise unsconstrained "wiggles" */
smv = xshapmag(p, v);
rv += smv;
/* Penalize if we have white > 1 or -ve black */
if (p->clipbw) {
double tp[3];
tp[0] = tp[1] = tp[2] = 1.0;
mxmfunc(p, v, xyz, tp);
if ((xyz[1] - 1.0) > err)
err = xyz[1] - 1.0;
tp[0] = tp[1] = tp[2] = 0.0;
mxmfunc(p, v, xyz, tp);
for (i = 0; i < 3; i++) {
if (-xyz[i] > err)
err = -v[i];
}
}
/* Penalize if we have -ve primaries */
if (p->clipprims) {
for (i = 0; i < 9; i++) {
if (-v[i] > err)
err = -v[i];
}
}
rv += err * 1000.0;
#ifdef DEBUG
if (g_deb)
printf("~9(%f)mxoptfunc returning %f\n",smv,rv);
#endif
return rv;
}
/* Matrix progress function handed to powell() */
static void mxprogfunc(void *pdata, int perc) {
mxopt *p = (mxopt *)pdata;
if (p->verb) {
printf("%c% 3d%%",cr_char,perc);
if (perc == 100)
printf("\n");
fflush(stdout);
}
}
/* Given a correction matrix, transform the matrix values */
static void mxtransform(mxopt *os, double mat[3][3]) {
double vec[3];
vec[0] = os->v[0]; vec[1] = os->v[3]; vec[2] = os->v[6];
icmMulBy3x3(vec, mat, vec);
os->v[0] = vec[0]; os->v[3] = vec[1]; os->v[6] = vec[2];
vec[0] = os->v[1]; vec[1] = os->v[4]; vec[2] = os->v[7];
icmMulBy3x3(vec, mat, vec);
os->v[1] = vec[0]; os->v[4] = vec[1]; os->v[7] = vec[2];
vec[0] = os->v[2]; vec[1] = os->v[5]; vec[2] = os->v[8];
icmMulBy3x3(vec, mat, vec);
os->v[2] = vec[0]; os->v[5] = vec[1]; os->v[8] = vec[2];
}
/* Setup and then return the optimized matrix fit in the mxopt structure. */
/* Return 0 on sucess, error code on failure. */
static int
createMatrix(
char *err, /* Return error message */
mxopt *os, /* Optimisation information */
int verb, /* NZ if verbose */
int nodp, /* Number of points */
cow *ipoints, /* Array of input points in XYZ space */
int isLab, /* nz if data points are Lab */
int quality, /* Quality metric, 0..3 (-1 == 2 orders only) */
int isLinear, /* NZ if pure linear, gamma = 1.0 */
int isGamma, /* NZ if gamma rather than shaper */
int isShTRC, /* NZ if shared TRCs */
int shape0gam, /* NZ if zero'th order shaper should be gamma function */
int clipbw, /* Prevent white > 1 and -ve black */
int clipprims, /* Prevent primaries going -ve */
double smooth, /* Smoothing factor (nominal 1.0) */
double scale /* Scale device values */
) {
double nweight = 1.0; /* Amount to weight neutral patches (make a parameter ?) */
int inputChan = 3; /* Must be RGB like */
int outputChan = 3; /* Must be the PCS */
int rsplflags = 0; /* Flags for scattered data rspl */
int e, f, i, j;
int maxits = 200; /* Optimisation stop params */
double stopon = 0.01; /* Absolute delta E change to stop on */
cow *points; /* Lab copy of ipoints */
int wix = -1; /* White patch index */
double rerr;
#ifdef DEBUG_PLOT
#define XRES 100
double xx[XRES];
double y1[XRES];
#endif /* DEBUG_PLOT */
if (verb)
rsplflags |= RSPL_VERBOSE;
/* Allocate the array passed to fit_rspl() */
if ((points = (cow *)malloc(sizeof(cow) * nodp)) == NULL) {
if (err != NULL)
sprintf(err,"Allocation of scattered coordinate array failed");
return 2;
}
/* Setup for optimising run */
if (verb != 0)
os->verb = verb;
else
os->verb = 0;
os->points = points;
os->nodp = nodp;
os->isShTRC = 0;
os->shape0gam = shape0gam;
os->smooth = smooth;
os->clipbw = clipbw;
os->clipprims = clipprims;
os->dscale = scale;
/* Set quality/effort factors */
if (quality >= 3) { /* Ultra high */
os->norders = 20;
maxits = 50000;
stopon = 1e-14;
} else if (quality == 2) { /* High */
os->norders = 12;
maxits = 5000;
stopon = 5e-6;
} else if (quality == 1) { /* Medium */
os->norders = 8;
maxits = 2000;
stopon = 5e-5;
} else if (quality == 0) { /* Low */
os->norders = 4;
maxits = 1000;
stopon = 5e-4;
} else { /* Ultra Low */
os->norders = 2;
maxits = 1000;
stopon = 5e-4;
}
if (os->norders > MXNORDERS)
os->norders = MXNORDERS;
/* Pick a white point for the real Lab conversion */
{
double wpy = -1e60;
double wp[3];
/* We assume that the input target is well behaved, */
/* and that it includes a white point patch, */
/* and that it has an extreme L value */
wix = -1;
for (i = 0; i < nodp; i++) {
double lab[3];
double yv;
/* Get the Lab value of this point */
if (isLab)
icmCpy3(lab, ipoints[i].v);
else
icmXYZ2Lab(&icmD50, lab, ipoints[i].v);
/* Tilt things towards D50 neutral white patches */
yv = lab[0] - 0.3 * sqrt(lab[1] * lab[1] + lab[2] * lab[2]);
if (yv > wpy) {
wpy = yv;
wix = i;
icmCpy3(wp, lab);
}
}
//printf("~1 picked point %d as Lab %f %f %f as white\n",wix,wp[0],wp[1],wp[2]);
icmLab2XYZ(&icmD50, wp, wp);
wp[0] /= wp[1]; /* Normalize */
wp[2] /= wp[1];
wp[1] = 1.0;
icmAry2XYZ(os->wp, wp);
/* We'll use this wp for delta E calculation when creating the matrix */
// if (os->verb) printf("Switching to L*a*b* white point %f %f %f\n",os->wp.X,os->wp.Y,os->wp.Z);
}
if (nweight < 1.0) /* Sanity */
nweight = 1.0;
/* Setup points ready for optimisation by converting */
/* data to Lab relative to our chosen white point */
for (i = 0; i < nodp; i++) {
double lch[3];
/* Copy input values */
for (e = 0; e < inputChan; e++)
points[i].p[e] = ipoints[i].p[e];
points[i].w = ipoints[i].w;
/* First copy and convert to XYZ */
if (isLab)
icmLab2XYZ(&icmD50, points[i].v, ipoints[i].v);
else
icmCpy3(points[i].v, ipoints[i].v);
/* Then to our wp relative Lab */
icmXYZ2Lab(&os->wp, points[i].v, points[i].v);
icmLab2LCh(lch, points[i].v);
/* Apply any neutral weighting */
if (lch[1] < 10.0) {
double w = nweight;
if (lch[1] > 5.0)
w = 1.0 + (nweight - 1.0) * (10.0 - lch[1])/(10.0 - 5.0);
points[i].w *= w;
}
//printf("~1 patch %d = Lab %f %f %f, C = %f w = %f\n",i,points[i].v[0], points[i].v[1], points[i].v[2], lch[1],points[i].w);
}
/* Set initial matrix optimisation values */
os->v[0] = 0.4; os->v[1] = 0.4; os->v[2] = 0.2; /* Matrix */
os->v[3] = 0.2; os->v[4] = 0.8; os->v[5] = 0.1;
os->v[6] = 0.02; os->v[7] = 0.15; os->v[8] = 1.3;
/* We try and take a homomorphic approach here, in an attempt */
/* to avoid getting trapped at a local minimum when a full */
/* set of shaper parameters are in play. */
/* Do a first pass just setting the matrix values */
os->isLinear = 1;
os->isGamma = 1;
os->optdim = 9;
os->v[9] = os->v[10] = os->v[11] = 1.0; /* Linear */
/* Set search area to starting values */
for (j = 0; j < os->optdim; j++)
os->sa[j] = 0.2; /* Matrix, Gamma, Offsets, harmonics */
if (os->verb)
printf("Creating matrix...\n");
if (powell(&rerr, os->optdim, os->v, os->sa, stopon, maxits,
mxoptfunc, (void *)os, mxprogfunc, (void *)os) != 0)
warning("Powell failed to converge, residual error = %f",rerr);
#ifndef NEVER
if (os->verb) {
printf("Matrix = %f %f %f\n",os->v[0], os->v[1], os->v[2]);
printf(" %f %f %f\n",os->v[3], os->v[4], os->v[5]);
printf(" %f %f %f\n",os->v[6], os->v[7], os->v[8]);
}
#endif /* NEVER */
/* Now optimize again with shaper or gamma curves */
if (!isLinear) {
double scgamma;
/* Start from linear, which is what was assumed for the matrix fit, */
/* and fit with a single shared gamma curve. */
os->isShTRC = 1;
os->isLinear = 0;
os->isGamma = 1;
os->optdim = 10;
os->v[9] = 1.0; /* Linear */
/* Set search area to starting values */
for (j = 0; j < os->optdim; j++)
os->sa[j] = 0.2; /* Matrix, Gamma, Offsets, harmonics */
if (os->verb)
printf("Creating matrix and single gamma curve...\n");
if (powell(&rerr, os->optdim, os->v, os->sa, stopon, maxits,
mxoptfunc, (void *)os, mxprogfunc, (void *)os) != 0)
warning("Powell failed to converge, residual error = %f",rerr);
scgamma = os->v[9];
if (isShTRC && !isGamma) {
#ifndef NEVER
if (os->verb) {
printf("Matrix = %f %f %f\n",os->v[0], os->v[1], os->v[2]);
printf(" %f %f %f\n",os->v[3], os->v[4], os->v[5]);
printf(" %f %f %f\n",os->v[6], os->v[7], os->v[8]);
printf("Gamma = %f\n",os->v[9]);
}
#endif /* NEVER */
/* Do final optimisation using full curve capability */
/* and fit first with a single shared curve. */
os->isShTRC = 1;
os->isLinear = 0;
os->isGamma = 0;
os->optdim = 9 + 2 + os->norders; /* Matrix, offset + orders */
os->v[9] = 0.0; /* Input offset */
os->v[10] = 0.0; /* Output offset */
if (shape0gam)
os->v[11] = 1.0; /* Gamma */
else
os->v[11] = 0.0; /* 0th Harmonic */
for (i = 12; i < os->optdim; i++)
os->v[i] = 0.0; /* Higher orders */
/* Set search area to starting values */
for (j = 0; j < os->optdim; j++)
os->sa[j] = 0.2; /* Matrix, Gamma, Offsets, harmonics */
if (os->verb)
printf("Creating matrix and single shaper curve...\n");
if (powell(&rerr, os->optdim, os->v, os->sa, stopon, maxits,
mxoptfunc, (void *)os, mxprogfunc, (void *)os) != 0)
warning("Powell failed to converge, residual error = %f",rerr);
scgamma = os->v[9];
}
/* For multiple curves, continue fitting */
if (!isShTRC) {
double mcgamma[3];
#ifndef NEVER
if (os->verb) {
printf("Matrix = %f %f %f\n",os->v[0], os->v[1], os->v[2]);
printf(" %f %f %f\n",os->v[3], os->v[4], os->v[5]);
printf(" %f %f %f\n",os->v[6], os->v[7], os->v[8]);
printf("Gamma = %f\n",os->v[9]);
}
#endif /* NEVER */
/* Fit matrix + multi gamma curves */
os->isShTRC = 0;
os->isLinear = 0;
os->isGamma = 1;
os->optdim = 12;
os->v[9] = os->v[10] = os->v[11] = scgamma; /* Single curve value */
/* Set search area to starting values */
for (j = 0; j < os->optdim; j++)
os->sa[j] = 0.2; /* Matrix, Gamma, Offsets, harmonics */
if (os->verb)
printf("Creating matrix and gamma curves...\n");
if (powell(&rerr, os->optdim, os->v, os->sa, stopon, maxits,
mxoptfunc, (void *)os, mxprogfunc, (void *)os) != 0)
warning("Powell failed to converge, residual error = %f",rerr);
mcgamma[0] = os->v[9];
mcgamma[1] = os->v[10];
mcgamma[2] = os->v[11];
if (!isGamma) {
#ifndef NEVER
if (os->verb) {
printf("Matrix = %f %f %f\n",os->v[0], os->v[1], os->v[2]);
printf(" %f %f %f\n",os->v[3], os->v[4], os->v[5]);
printf(" %f %f %f\n",os->v[6], os->v[7], os->v[8]);
printf("Gamma = %f %f %f\n",os->v[9], os->v[10], os->v[11]);
}
#endif /* NEVER */
/* Do final curves */
os->isShTRC = 0;
os->isLinear = 0;
os->isGamma = 0;
os->optdim = 9 + 6 + 3 * os->norders; /* Matrix, offset + orders */
os->v[9] = os->v[10] = os->v[11] = 0.0; /* Input offset */
os->v[12] = os->v[13] = os->v[14] = 0.0; /* Output offset */
if (shape0gam) {
os->v[15] = mcgamma[0];
os->v[16] = mcgamma[1];
os->v[17] = mcgamma[2];
} else
os->v[15] = os->v[16] = os->v[17] = 0.0; /* 0th Harmonic */
for (i = 18; i < os->optdim; i++)
os->v[i] = 0.0; /* Higher orders */
/* Set search area to starting values */
for (j = 0; j < os->optdim; j++)
os->sa[j] = 0.1; /* Matrix, Gamma, Offsets, harmonics */
if (os->verb)
printf("Creating matrix and curves...\n");
//g_deb = 1;
if (powell(&rerr, os->optdim, os->v, os->sa, stopon, maxits,
mxoptfunc, (void *)os, mxprogfunc, (void *)os) != 0)
warning("Powell failed to converge, residual error = %f",rerr);
}
}
}
if (os->clipprims) { /* Clip -ve primaries */
for (i = 0; i < 9; i++) {
if (os->v[i] < 0.0)
os->v[i] = 0.0;
}
}
#ifndef NEVER
if (os->verb) {
printf("Matrix = %f %f %f\n",os->v[0], os->v[1], os->v[2]);
printf(" %f %f %f\n",os->v[3], os->v[4], os->v[5]);
printf(" %f %f %f\n",os->v[6], os->v[7], os->v[8]);
if (!isLinear) { /* Creating input curves */
if (os->isGamma) { /* Creating input curves */
if (isShTRC)
printf("Gamma = %f\n",os->v[9]);
else
printf("Gamma = %f %f %f\n",os->v[9], os->v[10], os->v[11]);
} else { /* Creating input curves */
if (isShTRC) {
printf("Input offset = %f\n",os->v[9]);
printf("Output offset = %f\n",os->v[10]);
} else {
printf("Input offset = %f %f %f\n",os->v[9], os->v[10], os->v[11]);
printf("Output offset = %f %f %f\n",os->v[12], os->v[13], os->v[14]);
}
for (j = 0; j < os->norders; j++) {
if (isShTRC) {
if (shape0gam && j == 0)
printf("gamma = %f\n", os->v[11 + j]);
else
printf("%d harmonics = %f\n",j, os->v[11 + j]);
} else {
if (shape0gam && j == 0)
printf("%d gamma = %f %f %f\n",j, os->v[15 + j * 3],
os->v[16 + j * 3], os->v[17 + j * 3]);
else
printf("%d harmonics = %f %f %f\n",j, os->v[15 + j * 3],
os->v[16 + j * 3], os->v[17 + j * 3]);
}
}
}
}
}
#endif /* NEVER */
#ifdef NEVER /* Check DE of fit */
{
double xyz[3], txyz[3];
for (i = 0; i < nodp; i++) {
mxmfunc(os, os->v, xyz, ipoints[i].p);
if (isLab)
icmLab2XYZ(&icmD50, txyz, ipoints[i].v);
else
icmCpy3(txyz, ipoints[i].v);
printf("~1 point %d DE %f\n", i, icmXYZLabDE(&icmD50, txyz, xyz));
}
}
#endif
/* Free the coordinate lists */
free(points);
return 0;
}
/* Apply a chromatic transform to the matrix to force the given */
/* xyz value (typically white) to be exact */
static void icxMM_force_exact(icxMatrixModel *p, double *targ, double *rgb) {
mxopt *os = (mxopt *)p->imp;
double txyz[3], axyz[3]; /* Target & actual xyz */
icmXYZNumber _tp, _ap;
double cmat[3][3]; /* Model transform matrix */
if (p->isLab)
icmLab2XYZ(&icmD50, txyz, targ);
else
icmCpy3(txyz, targ);
mxmfunc(os, os->v, axyz, rgb);
icmAry2XYZ(_ap, axyz);
icmAry2XYZ(_tp, txyz);
if (p->picc != NULL)
p->picc->chromAdaptMatrix(p->picc, ICM_CAM_NONE, NULL, cmat, _tp, _ap);
else
icmChromAdaptMatrix(ICM_CAM_BRADFORD, _tp, _ap, cmat);
/* Apply correction to fine tune matrix. */
mxtransform(os, cmat);
}
static void icxMM_lookup(icxMatrixModel *p, double *out, double *in) {
mxopt *os = (mxopt *)p->imp;
mxmfunc(os, os->v, out, in);
if (p->isLab)
icmXYZ2Lab(&icmD50, out, out);
}
static void icxMM_del(icxMatrixModel *p) {
free(p->imp);
free(p);
}
/* Create a matrix model of a set of points, and return an object to lookup */
/* points from the model. Return NULL on error. */
icxMatrixModel *new_MatrixModel(
icc *picc, /* ICC profile used to set cone space matrix, NULL for Bradford. */
int verb, /* NZ if verbose */
int nodp, /* Number of points */
cow *ipoints, /* Array of input points in XYZ space */
int isLab, /* nz if data points are Lab */
int quality, /* Quality metric, 0..3 (-1 == 2 orders only) */
int isLinear, /* NZ if pure linear, gamma = 1.0 */
int isGamma, /* NZ if gamma rather than shaper */
int isShTRC, /* NZ if shared TRCs */
int shape0gam, /* NZ if zero'th order shaper should be gamma function */
int clipbw, /* Prevent white > 1 and -ve black */
int clipprims, /* Prevent primaries going -ve */
double smooth, /* Smoothing factor (nominal 1.0) */
double scale /* Scale device values */
) {
icxMatrixModel *p;
if ((p = (icxMatrixModel *) calloc(1,sizeof(icxMatrixModel))) == NULL)
return NULL;
p->picc = picc;
p->force = icxMM_force_exact;
p->lookup = icxMM_lookup;
p->del = icxMM_del;
if ((p->imp = (void *) calloc(1,sizeof(mxopt))) == NULL) {
free(p);
return NULL;
}
if (createMatrix(NULL, (mxopt *)p->imp, verb, nodp, ipoints, isLab, quality,
isLinear, isGamma, isShTRC, shape0gam,
clipbw, clipprims, smooth, scale) != 0) {
free(p->imp);
free(p);
return NULL;
}
p->isLab = isLab;
return p;
}
/* Callback function for create_matrix_xforms() */
static void matrix_xforms_infunc(void *cbctx, double *out, double *in, int tn) {
mxopt *p = (mxopt *)cbctx;
int j;
for (j = 0; j < 3; j++)
mxmfunc1(p, j, p->v, out + j, in + j);
}
/* Create icxLuMatrix and undelying tone reproduction curves and */
/* colorant tags from the supplied scattered data points. */
/* The scattered data is assumed to map Device -> native PCS (ie. dir = Fwd) */
/* Set error code and message in icco on error. */
/* Return the error code */
int set_icxLuMatrix(
icc *icco,
int flags, /* white/black point flags */
icxTransformCreateType icctype, /* Profile creation type */
int nodp, /* Number of points */
int nodpbw, /* Number of points to look for white & black patches in */
cow *ipoints, /* Array of input points in XYZ space */
icxMatrixModel *skm, /* Optional skeleton model (not used here) */
double dispLuminance, /* > 0.0 if display luminance value and is known */
double wpscale, /* > 0.0 if input white point is to be scaled */
int quality, /* Quality metric, 0..3 */
double smooth, /* Curve smoothing, nominally 1.0 */
int isShTRC, /* Matrix - NZ if shared TRCs */
int isGamma, /* Matrix - NZ if gamma rather than shaper */
int isLinear, /* Matrix - NZ if pure linear, gamma = 1.0 */
int inputEnt, /* Num of in-table entries */
int inv_inputEnt /* Num of inverse in-table entries */
) {
icmHeader *h = icco->header; /* Pointer to icc header */
int inputChan = 3; /* Must be RGB like */
int outputChan = 3; /* Must be the PCS */
int rsplflags = 0; /* Flags for scattered data rspl */
int e, f, i, j;
int maxits = 200; /* Optimisation stop params */
double stopon = 0.01; /* Absolute delta E change to stop on */
mxopt os; /* Optimisation information */
/* If ICX_SET_WHITE | ICX_SET_BLACK: */
double wp[3]; /* Absolute White point in XYZ */
double bp[3]; /* Absolute Black point in XYZ */
double dw[MXDI]; /* Device white value to adjust to be D50 */
double db[MXDI]; /* Device balck value */
double dgw[3]; /* Device space gamut boundary white for ICX_SET_WHITE_US */
double fromAbs[3][3]; /* From abs to relative */
double toAbs[3][3]; /* To abs from relative */
cow *rpoints = NULL; /* Aprox. relative in->output values */
#ifdef DEBUG_PLOT
#define XRES 100
double xx[XRES];
double y1[XRES];
#endif /* DEBUG_PLOT */
if (flags & ICX_VERBOSE)
rsplflags |= RSPL_VERBOSE;
/* Sanity check for matrix */
if (h->pcs != icSigXYZData) {
return icm_err(icco, ICX_ERR_BAD_PCS, "Can't create matrix profile with PCS of %s !",
icm2str(icmColorSpaceSig, h->pcs));
}
/* ------------------------------- */
/* Choose a white and black point */
if (flags & (ICX_SET_WHITE | ICX_SET_BLACK)) {
if (flags & ICX_VERBOSE)
printf("Find white & black points\n");
/* Compute device white and black points as if */
/* we are doing an Output or Display device */
{
switch (h->colorSpace) {
case icSigCmyData:
for (e = 0; e < inputChan; e++) {
dw[e] = 0.0;
db[e] = 1.0;
}
break;
case icSigRgbData:
for (e = 0; e < inputChan; e++) {
dw[e] = 1.0;
db[e] = 0.0;
}
break;
default:
return icm_err(icco,ICX_ERR_BAD_DEV,
"set_icxLuMatrix: can't handle color space %s",
icm2str(icmColorSpaceSig, h->colorSpace));
}
}
/* dw is what we want for dgw[], used for XFIT_OUT_WP_REL_US */
for (e = 0; e < inputChan; e++)
dgw[e] = dw[e];
/* If this is actuall an input device, lookup wp & bp */
/* and override dwhite & dblack */
if (h->deviceClass == icSigInputClass) {
double wpy = -1e60, bpy = 1e60;
int wix = -1, bix = -1;
/* We assume that the input target is well behaved, */
/* and that it includes a white and black point patch, */
/* and that they have the extreme L/Y values */
/*
NOTE that this may not be the best approach !
It may be better to average the chromaticity
of all the neutral seeming patches, since
the whitest patch may have (for instance)
a blue tint.
*/
/* Discover the white and black patches */
for (i = 0; i < nodpbw; i++) {
double labv[3], yv;
/* Create D50 Lab to allow some chromatic sensitivity */
/* in picking the white point */
icmXYZ2Lab(&icmD50, labv, ipoints[i].v);
#ifdef NEVER
/* Choose Y */
if (ipoints[i].v[1] > wpy) {
wp[0] = ipoints[i].v[0];
wp[1] = ipoints[i].v[1];
wp[2] = ipoints[i].v[2];
for (e = 0; e < inputChan; e++)
dw[e] = ipoints[i].p[e];
wpy = ipoints[i].v[1];
wix = i;
}
#else
/* Tilt things towards D50 neutral white patches */
yv = labv[0] - 0.3 * sqrt(labv[1] * labv[1] + labv[2] * labv[2]);
if (yv > wpy) {
wp[0] = ipoints[i].v[0];
wp[1] = ipoints[i].v[1];
wp[2] = ipoints[i].v[2];
for (e = 0; e < inputChan; e++)
dw[e] = ipoints[i].p[e];
wpy = yv;
wix = i;
}
#endif
if (ipoints[i].v[1] < bpy) {
bp[0] = ipoints[i].v[0];
bp[1] = ipoints[i].v[1];
bp[2] = ipoints[i].v[2];
for (e = 0; e < inputChan; e++)
db[e] = ipoints[i].p[e];
bpy = ipoints[i].v[1];
bix = i;
}
}
if (flags & ICX_VERBOSE) {
printf("Picked white patch %d with dev = %s\n XYZ = %s, Lab = %s\n",
wix+1, icmPdv(inputChan, dw), icmPdv(3, wp), icmPLab(wp));
printf("Picked black patch %d with dev = %s\n XYZ = %s, Lab = %s\n",
bix+1, icmPdv(inputChan, db), icmPdv(3, bp), icmPLab(bp));
}
} else {
/* We assume that the display target is well behaved, */
/* and that it includes a white point patch. */
int nw = 0;
wp[0] = wp[1] = wp[2] = 0.0;
switch (h->colorSpace) {
case icSigCmyData:
for (i = 0; i < nodpbw; i++) {
if (ipoints[i].p[0] < 0.001
&& ipoints[i].p[1] < 0.001
&& ipoints[i].p[2] < 0.001) {
wp[0] += ipoints[i].v[0];
wp[1] += ipoints[i].v[1];
wp[2] += ipoints[i].v[2];
nw++;
}
}
break;
case icSigRgbData:
for (i = 0; i < nodpbw; i++) {
if (ipoints[i].p[0] > 0.999
&& ipoints[i].p[1] > 0.999
&& ipoints[i].p[2] > 0.999) {
wp[0] += ipoints[i].v[0];
wp[1] += ipoints[i].v[1];
wp[2] += ipoints[i].v[2];
nw++;
}
}
break;
default:
return icm_err(icco,ICX_ERR_BAD_DEV,
"set_icxLuMatrix: can't handle color space %s",
icm2str(icmColorSpaceSig, h->colorSpace));
}
if (nw == 0) {
return icm_err(icco, ICX_ERR_NO_WP,
"set_icxLuMatrix: can't handle test points without a white patch");
}
wp[0] /= (double)nw;
wp[1] /= (double)nw;
wp[2] /= (double)nw;
if (flags & ICX_VERBOSE) {
printf("Initial white point = %f %f %f\n",wp[0],wp[1],wp[2]);
}
/* Need to lookup bp[] before we set the tag */
}
/* Create some abs<->rel chromatic conversions */
{
icmXYZNumber _wp;
icmAry2XYZ(_wp, wp);
/* Absolute->Aprox. Relative Adaptation matrix, and */
/* Aproximate relative to absolute conversion matrix */
icco->chromAdaptMatrix(icco, ICM_CAM_NONE, toAbs, fromAbs, icmD50, _wp);
}
} else {
icmSetUnity3x3(fromAbs);
icmSetUnity3x3(toAbs);
icmCpy3(wp, icmD50_ary3);
}
/* Force do nothing/Absolute white point */
if ((flags & ICX_SET_WHITE_ABS) == ICX_SET_WHITE_ABS) {
icmSetUnity3x3(fromAbs);
icmSetUnity3x3(toAbs);
icmCpy3(wp, icmD50_ary3);
}
/* Create copy of input points with output converted to white relative */
/* Allow one extra point for possible bpo value */
if ((rpoints = (cow *)malloc((nodp+1) * sizeof(cow))) == NULL)
return icm_err(icco, ICX_ERR_MALLOC, "set_icxLuMatrix: malloc failed");
for (i = 0; i < nodp; i++) {
rpoints[i].w = ipoints[i].w;
for (e = 0; e < inputChan; e++)
rpoints[i].p[e] = ipoints[i].p[e];
for (f = 0; f < outputChan; f++)
rpoints[i].v[f] = ipoints[i].v[f];
/* abs out -> aprox. rel out */
icmMulBy3x3(rpoints[i].v, fromAbs, rpoints[i].v);
}
/* ------------------------------- */
/* (Use a gamma curve as 0th order shape) */
if ((icco->e.c = createMatrix(icco->e.m, &os, flags & ICX_VERBOSE ? 1 : 0,
nodp, rpoints, 0, quality,
isLinear, isGamma, isShTRC, 1,
flags & ICX_CLIP_WB ? 1 : 0,
flags & ICX_CLIP_PRIMS ? 1 : 0,
smooth, 1.0)) != ICX_ERR_OK) {
free(rpoints);
return icco->e.c;
}
free(rpoints); rpoints = NULL;
/* 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 ((flags & ICX_SET_WHITE) != 0
&& (flags & ICX_SET_WHITE_ABS) != ICX_SET_WHITE_ABS) {
double aw[3]; /* aprox rel. white */
icmXYZNumber _wp; /* Uncorrected dw maps to _wp */
double cmat[3][3]; /* Model correction matrix */
if (flags & ICX_VERBOSE)
printf("Doing White point fine tune:\n");
/* See what the aprox. relative white point has turned out to be, */
/* by looking up the device white in the current conversion */
mxmfunc(&os, os.v, aw, dw);
if (flags & ICX_VERBOSE) {
printf("Before fine tune, rel WP = XYZ %s, Lab %s\n", icmPdv(3,aw), icmPLab(aw));
}
/* Matrix needed to correct aprox white to target D50 */
icmAry2XYZ(_wp, aw); /* Aprox relative target white point */
icco->chromAdaptMatrix(icco, ICM_CAM_NONE, NULL, cmat, icmD50, _wp); /* Correction */
/* Compute the current absolute white point */
icmMulBy3x3(wp, toAbs, aw);
/* Apply correction to fine tune matrix. */
mxtransform(&os, cmat);
/* Fix relative conversions to leave absolute response unchanged. */
icmAry2XYZ(_wp, wp); /* Actual white point */
icco->chromAdaptMatrix(icco, ICM_CAM_NONE, toAbs, fromAbs, icmD50, _wp);
if (flags & ICX_VERBOSE) {
double tw[3];
mxmfunc(&os, os.v, tw, dw); /* Lookup white again */
printf("After fine tune, rel WP = XYZ %s, Lab %s\n", icmPdv(3, tw), icmPLab(tw));
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 (flags & ICX_VERBOSE) {
printf("White manual point scale %f\n", wpscale);
}
}
/* If we are going to auto scale the WP to avoid clipping */
/* values above the WP: (not so important for matrix profiles ?) */
if ((flags & ICX_SET_WHITE_US) == ICX_SET_WHITE_US) {
double tw[3], 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 */
mxmfunc(&os, os.v, tw, dgw);
icmMulBy3x3(tw, toAbs, tw); /* Convert to absolute */
mxY = tw[1];
icmCpy3(bw, tw);
//printf("~1 1,1,1 Y = %f\n",tw[1]);
/* See what the device white point value scaled to 1 produces */
mxd = -1.0;
for (e = 0; e < inputChan; e++) {
if (dw[e] > mxd)
mxd = dw[e];
}
for (e = 0; e < inputChan; e++)
ndw[e] = dw[e]/mxd;
mxmfunc(&os, os.v, tw, ndw);
icmMulBy3x3(tw, toAbs, tw); /* Convert to absolute */
//printf("~1 ndw = %f %f %f Y = %f\n",ndw[0],ndw[1],ndw[2],tw[1]);
if (tw[1] > mxY) {
mxY = tw[1];
icmCpy3(bw, tw);
}
/* Compute WP scale factor needed to fit mxY */
if (mxY > wp[1]) {
uswpscale = mxY/wp[1];
wpscale *= uswpscale;
if (flags & ICX_VERBOSE) {
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 (flags & ICX_CLIP_WB) {
if ((wp[1] * wpscale) > 1.0) {
wpscale = 1.0/wp[1]; /* Make wp Y = 1.0 */
if (flags & ICX_VERBOSE) {
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;
double cmat[3][3]; /* Model correction matrix */
/* Create inverse scaling matrix for relative rspl data */
icmSetUnity3x3(cmat);
icmScale3x3(cmat, cmat, 1.0/wpscale);
/* Inverse scale the matrix */
mxtransform(&os, cmat);
/* Scale the WP */
icmScale3(wp, wp, wpscale);
/* Fix absolute conversions to leave absolute response unchanged. */
icmAry2XYZ(_wp, wp); /* Actual white point */
icco->chromAdaptMatrix(icco, ICM_CAM_NONE, toAbs, fromAbs, icmD50, _wp);
}
/* Look up the actual black point */
if (flags & ICX_SET_BLACK) {
/* Look black point up in dev->rel model */
mxmfunc(&os, os.v, bp, db);
/* Convert from relative to Absolute colorimetric */
icmMulBy3x3(bp, toAbs, bp);
/* Got XYZ black point in bp[] */
if (flags & ICX_VERBOSE) {
printf("Black point XYZ = %s, Lab = %s\n", icmPdv(3,bp),icmPLab(bp));
}
if (flags & ICX_CLIP_WB) {
if (bp[0] < 0.0 || bp[1] < 0.0 || bp[1] < 0.0) {
if (bp[0] < 0.0)
bp[0] = 0.0;
if (bp[1] < 0.0)
bp[1] = 0.0;
if (bp[2] < 0.0)
bp[2] = 0.0;
if (flags & ICX_VERBOSE)
printf("Black point clipped to XYZ = %s, Lab = %s\n",icmPdv(3,bp),icmPLab(bp));
}
}
}
if (flags & (ICX_SET_WHITE | ICX_SET_BLACK)) {
/* If this is a display, adjust the absolute white point to be */
/* exactly Y = 1.0, and compensate the matrix, dispLuminance */
/* and black point accordingly. */
if (h->deviceClass == icSigDisplayClass) {
double cmat[3][3]; /* Model correction matrix */
double scale = 1.0/wp[1];
if (flags & ICX_VERBOSE)
printf("Scaling White Point by %f to make Y = 1.0\n", scale);
/* Scale the WP & BP*/
icmScale3(wp, wp, scale);
icmScale3(bp, bp, scale);
/* Inverse scale the luminance */
dispLuminance /= scale;
}
/* Absolute luminance tag */
if (flags & ICX_WRITE_WBL
&& h->deviceClass == icSigDisplayClass
&& dispLuminance > 0.0) {
icmXYZArray *wo;
/* Delete any existing tag */
if (icco->delete_tag_quiet(icco, icSigLuminanceTag) != ICM_ERR_OK)
return icco->e.c;
/* Note that tag types icSigXYZType and icSigXYZArrayType are identical */
if ((wo = (icmXYZArray *)icco->add_tag(
icco, icSigLuminanceTag, icSigXYZArrayType)) == NULL)
return icco->e.c;
wo->count = 1;
wo->allocate(wo); /* Allocate space */
wo->data[0].X = 0.0;
wo->data[0].Y = dispLuminance;
wo->data[0].Z = 0.0;
if (flags & ICX_VERBOSE)
printf("Display Luminance = %f\n", wo->data[0].Y);
}
/* Add the white and black tags */
if ((flags & ICX_WRITE_WBL)
&& (flags & ICX_SET_WHITE)) { /* White Point Tag: */
icmXYZArray *wo;
/* Delete any existing tag */
if (icco->delete_tag_quiet(icco, icSigMediaWhitePointTag) != ICM_ERR_OK)
return icco->e.c;
/* Note that tag types icSigXYZType and icSigXYZArrayType are identical */
if ((wo = (icmXYZArray *)icco->add_tag(
icco, icSigMediaWhitePointTag, icSigXYZArrayType)) == NULL)
return icco->e.c;
wo->count = 1;
wo->allocate(wo); /* Allocate space */
wo->data[0].X = wp[0];
wo->data[0].Y = wp[1];
wo->data[0].Z = wp[2];
if (flags & ICX_VERBOSE)
printf("White point XYZ = %f %f %f\n",wp[0],wp[1],wp[2]);
}
if ((flags & ICX_WRITE_WBL)
&& (flags & ICX_SET_BLACK)) { /* Black Point Tag: */
icmXYZArray *wo;
/* Delete any existing tag */
if (icco->delete_tag_quiet(icco, icSigMediaBlackPointTag) != ICM_ERR_OK)
return icco->e.c;
/* Note that tag types icSigXYZType and icSigXYZArrayType are identical */
if ((wo = (icmXYZArray *)icco->add_tag(
icco, icSigMediaBlackPointTag, icSigXYZArrayType)) == NULL)
return icco->e.c;
wo->count = 1;
wo->allocate(wo); /* Allocate space */
wo->data[0].X = bp[0];
wo->data[0].Y = bp[1];
wo->data[0].Z = bp[2];
if (flags & ICX_VERBOSE)
printf("Black point XYZ = %f %f %f\n",bp[0],bp[1],bp[2]);
}
if (flags & ICX_CLIP_PRIMS) {
for (i = 0; i < 9; i++) {
if (os.v[i] < 0.0)
os.v[i] = 0.0;
}
}
}
if (flags & ICX_VERBOSE)
printf("Done gamma/shaper and matrix creation\n");
{
int nsigs = 1;
icmXformSigs sigs[2];
double mx[3][3]; /* [outn][inn] Matrix values */
double *gamma = NULL; /* Gamma values if isGamma */
switch (icctype) {
default:
sigs[0].sig = icmSigShaperMatrix;
sigs[0].ttype = icmSigShaperMatrixType;
/* Default ICC Version used */
break;
}
for (j = e = 0; j < 3; j++) {
for (i = 0; i < 3; i++, e++) {
mx[j][i] = os.v[e];
}
}
if (isGamma)
gamma = &os.v[9]; /* Gamma values */
/* Use icclib helper function to create tags and set their values */
/* Returns ec */
if (icco->create_matrix_xforms(
icco, ICM_CREATE_FLAG_NONE,
(void *)&os, /* Opaque callback context pointer value */
nsigs, sigs, /* Signatures and tag types for each table */
(unsigned int)inputEnt, (unsigned int)inv_inputEnt, /* fwd & rev Table resolutions */
h->colorSpace, h->pcs, /* ina & out color spaces */
matrix_xforms_infunc, /* Callback function */
mx, /* [outn][inn] Matrix values */
NULL, /* Constants if any */
isShTRC, gamma, isLinear /* curve flags */
) != 0) {
if (flags & ICX_VERBOSE)
printf("Matrix profile creation failed\n");
return icco->e.c;
}
}
if (flags & ICX_VERBOSE)
printf("Profile done\n");
return ICX_ERR_OK;
}
/* ========================================================= */
/* Given an xicc lookup object, returm a gamut object. */
/* Note that the PCS must be Lab or Jab */
/* Return NULL on error, check e.c+err for reason */
static gamut *icxLuMatrixGamut(
icxLuBase *plu, /* this */
double detail /* gamut detail level, 0.0 = def */
) {
xicc *p = plu->pp; /* parent xicc */
icxLuMatrix *lumat = (icxLuMatrix *)plu; /* Lookup xMatrix type object */
icColorSpaceSignature pcs;
icmLookupFunc func;
double white[3], black[3], kblack[3];
gamut *gam;
int res; /* Sample point resolution */
int i, e;
if (detail == 0.0)
detail = 10.0;
/* get some details */
plu->spaces(plu, NULL, NULL, NULL, NULL, NULL, NULL, &func, &pcs);
if (func != icmFwd && func != icmBwd) {
p->e.c = 1;
sprintf(p->e.m,"Creating Gamut surface for anything other than Device <-> PCS is not supported.");
return NULL;
}
if (pcs != icSigLabData && pcs != icxSigJabData) {
p->e.c = 1;
sprintf(p->e.m,"Creating Gamut surface PCS of other than Lab or Jab is not supported.");
return NULL;
}
gam = new_gamut(detail, pcs == icxSigJabData, 0);
/* Explore the gamut by itterating through */
/* it with sample points in device space. */
res = (int)(600.0/detail); /* Establish an appropriate sampling density */
if (res < 40)
res = 40;
/* Since matrix profiles can't be non-monotonic, */
/* just itterate through the surface colors. */
for (i = 0; i < 3; i++) {
int co[3];
int ep[3];
int co_e = 0;
for (e = 0; e < 3; e++) {
co[e] = 0;
ep[e] = res;
}
ep[i] = 2;
while (co_e < 3) {
double in[3];
double out[3];
for (e = 0; e < 3; e++) /* Convert count to input value */
in[e] = co[e]/(ep[e]-1.0);
/* Always use the device->PCS conversion */
if (lumat->fwd_lookup((icxLuBase *)lumat, out, in) > 1)
error ("%d, %s",p->e.c,p->e.m);
gam->expand(gam, out);
/* Increment the counter */
for (co_e = 0; co_e < 3; co_e++) {
co[co_e]++;
if (co[co_e] < ep[co_e])
break; /* No carry */
co[co_e] = 0;
}
}
}
#ifdef NEVER
/* Try it twice */
for (i = 0; i < 3; i++) {
int co[3];
int ep[3];
int co_e = 0;
for (e = 0; e < 3; e++) {
co[e] = 0;
ep[e] = res;
}
ep[i] = 2;
while (co_e < 3) {
double in[3];
double out[3];
for (e = 0; e < 3; e++) /* Convert count to input value */
in[e] = co[e]/(ep[e]-1.0);
/* Always use the device->PCS conversion */
if (lumat->fwd_lookup((icxLuBase *)lumat, out, in) > 1)
error ("%d, %s",p->e.c,p->err);
gam->expand(gam, out);
/* Increment the counter */
for (co_e = 0; co_e < 3; co_e++) {
co[co_e]++;
if (co[co_e] < ep[co_e])
break; /* No carry */
co[co_e] = 0;
}
}
}
#endif
#ifdef NEVER // (doesn't seem to make much difference)
/* run along the primary ridges in more detail too */
/* just itterate through the surface colors. */
for (i = 0; i < 3; i++) {
int j;
double in[3];
double out[3];
res *= 4;
for (j = 0; j < res; j++) {
double vv = i/(res-1.0);
in[0] = in[1] = in[2] = vv;
in[i] = 0.0;
if (lumat->fwd_lookup((icxLuBase *)lumat, out, in) > 1)
error ("%d, %s",p->e.c,p->err);
gam->expand(gam, out);
in[0] = in[1] = in[2] = 0.0;
in[i] = vv;
if (lumat->fwd_lookup((icxLuBase *)lumat, out, in) > 1)
error ("%d, %s",p->e.c,p->err);
gam->expand(gam, out);
}
}
#endif
/* Put the white and black points in the gamut */
plu->efv_wh_bk_points(plu, white, black, kblack);
gam->setwb(gam, white, black, kblack);
/* set the cusp points by itterating through the 0 & 100% colorant combinations */
{
DCOUNT(co, 3, 3, 0, 0, 2);
gam->setcusps(gam, 0, NULL);
DC_INIT(co);
while(!DC_DONE(co)) {
int e;
double in[3];
double out[3];
if (!(co[0] == 0 && co[1] == 0 && co[2] == 0)
&& !(co[0] == 1 && co[1] == 1 && co[2] == 1)) { /* Skip white and black */
for (e = 0; e < 3; e++)
in[e] = (double)co[e];
/* Always use the device->PCS conversion */
if (lumat->fwd_lookup((icxLuBase *)lumat, out, in) > 1)
error ("%d, %s",p->e.c,p->e.m);
gam->setcusps(gam, 3, out);
}
DC_INC(co);
}
gam->setcusps(gam, 2, NULL);
}
#ifdef NEVER /* Not sure if this is a good idea ?? */
gam->getwb(gam, NULL, NULL, white, black); /* Get the actual gamut white and black points */
gam->setwb(gam, white, black); /* Put it back as colorspace one */
#endif
return gam;
}
#ifdef DEBUG
#undef DEBUG
#endif
Reference in New Issue View Git Blame Copy Permalink