diff options
Diffstat (limited to 'xicc/xlut.c')
| -rw-r--r-- | xicc/xlut.c | 4271 |
1 files changed, 4271 insertions, 0 deletions
diff --git a/xicc/xlut.c b/xicc/xlut.c new file mode 100644 index 0000000..21c286e --- /dev/null +++ b/xicc/xlut.c @@ -0,0 +1,4271 @@ + +/* + * International Color Consortium color transform expanded support + * + * Author: Graeme W. Gill + * Date: 2/7/00 + * Version: 1.00 + * + * Copyright 2000, 2001 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. + * + * This is the third major version of xlut.c (originally called xlut2.c) + * Based on the old xlut.c code (preserved as xlut1.c) + * This version uses xfit.c to do the curve and rspl fitting. + */ + +/* + * This module provides the expanded icclib functionality + * for lut based profiles. + * This file is #included in xicc.c, to keep its functions private. + * + * This version creates both input and output 1D luts by + * optimising the accuracy of the profile for a linear clut. + * + */ + +/* + TTBD: + + See gamut/gammap.c for notes on the desirability of + a non-minimum delta E colorimetric intent as the default. + + Should fix gamut's so that they have enough information + to spefify a gamut mapping without the need for + a source colorspace profile, and fix the code + to work with just a gamut. This would take us further + towards supporting the PRMG reference gamut interoperability + as an option. + + Should the input profile white point determination + be made a bit smarter about determining the chromaticity ? + ie. not take on the tint of the whitest patch, but an + average of neutral patches ? + + Should xlutfix.c be revived (also adding ICM_CLUT_SET_APXLS support), + to improve "bumpy black" problem ? + */ + +/* + + NOTE :- an alternative to the way display profile absolute is handled here + would be to always chromatically adapt the illuminant to D50, and encode + that in the Chromatic adapation tag. To make absolute colorimetric + do anything useful though, the chromatic adapation tag would + have to be used for absolute intent. + This may be the way of improving compatibility with other systems, + but would break compatibility with existing Argyll profiles, + unless special measures are taken: + + ie. + + 1) if (display profile & using chromatic adaptation tag) + Create Bradford chromatic adapation matrix and store it in tag + Adapt all the readings using Bradford + Create white point and store it in tag (white point will be D50) + Adapt all the readings to the white point using wrong Von-Kries (no change) + Store relative colorimetric cLUT + Set version >= 2.4 + + else + 2) if (display scheme A or using Argyll historical printer scheme) + Create white point and store it in tag + Adapt all the readings to the white point using Bradford + Store relative colorimetric tag + Set version < 2.4 for V2 profile + Add private Absolute Transform Matrix (labels V4 profile) + + 3) else (display scheme B or strict ICC printer compatibility) + Create white point and store it in tag + Adapt all the readings to the white point using Wrong Von-Kries + Store relative colorimetric tag + Set version >= 2.4 + + + Argyll Processing for each type + + 1) if display and chromatic adapation matrix + Un-adapt matrix or cLUT using wrong Von-Kries from white point + Un-adapt matrix or cLUT using chromatic matrix + Un-adapt apparant white point & black point using chromatic transform + + if (not absolute intent) + Create Bradford transfor from white to PCS D50 + Adapt all matrix or cLUT + else + 2) if (display scheme A or using Argyll < V2.4. profile + or find Absolute Transform Matrix) + if (absolute intent) + Un-adapt matrix or cLUT using Bradford from white point + + 3) else (display scheme B or !Argyll profile or ( >= V2.4 profile + and !Absolute Transform Matrix)) + Un-adapt matrix or cLUT using wrong Von-Kries from white point + + if (not absolute intent) + Create Bradford transfor from white to PCS D50 + Adapt all matrix or cLUT to white + + + The problem with this is that it wouldn't do the right thing on old Argyll + type profiles that weren't labeled or recognized. + + Is there a way of recognizing Bradford Absolute transform Matricies if + the color chromaticities are given ? + + */ + +/* + A similar condrum is that it seems that an unwritten convention for + V2 profiles is to scale the black point of the perceptual and + saturation tables to 0 (Part of the V4 spec is to scale to Y = 3.1373). + + To get better gamut mapping we should therefore unscale the perceptual + and saturation A2B table to have the same black point as the colorimetric + table before computing the gamut mapping, and then apply the opposite + transform to our perceptual B2A and A2B tables. + + */ + +/* + There is interesting behaviour for Jab 0,0,0 in, in that + it gets mapped to (something like) Lab -1899.019855 -213.574625 -231.914285 + which after per-component clipping of the inv out tables looks like + 0, -128, -128, which may be mapped to (say) RGB 0.085455 1.000000 0.936951, + which is not black. + + */ + +#include "xfit.h" + +#undef USE_CIE94_DE /* [Undef] Use CIE94 delta E measure when creating in/out curves */ + /* Don't use CIE94 because it makes peak error worse ? */ + +#undef DEBUG /* [Undef] Verbose debug information */ +#undef CNDTRACE /* [Undef] Enable xicc->trace conditional verbosity */ +#undef DEBUG_OLUT /* [Undef] Print steps in overall fwd & rev lookups */ +#undef DEBUG_RLUT /* [Undef] Print values being reverse lookup up */ +#undef DEBUG_SPEC /* [Undef] Debug some specific cases */ +#undef INK_LIMIT_TEST /* [Undef] Turn off input tables for ink limit testing */ +#undef CHECK_ILIMIT /* [Undef] Do sanity checks on meeting ink limit */ +#undef WARN_CLUT_CLIPPING /* [Undef] Print warning if setting clut clips */ +#undef DISABLE_KCURVE_FILTER /* [Undef] don't filter the Kcurve */ +#undef REPORT_LOCUS_SEGMENTS /* [Undef[ Examine how many segments there are in aux inversion */ + +#define SHP_SMOOTH 1.0 /* Input shaper curve smoothing */ +#define OUT_SMOOTH1 1.0 /* Output shaper curve smoothing for L*, X,Y,Z */ +#define OUT_SMOOTH2 1.0 /* Output shaper curve smoothing for a*, b* */ + +#define CAMCLIPTRANS 1.0 /* [1.0] Cam clipping transition region Delta E */ +#define USECAMCLIPSPLINE /* [def] use spline blend */ + +/* + * TTBD: + * + * Reverse lookup of Lab + * Make NEARCLIP the default ?? + * + * XYZ vector clipping isn't implemented. + * + * Some of the error handling is crude. Shouldn't use + * error(), should return status. + */ + +#ifndef _CAT2 +#define _CAT2(n1,n2) n1 ## n2 +#define CAT2(n1,n2) _CAT2(n1,n2) +#endif + + + +static double icxLimitD(icxLuLut *p, double *in); /* For input' */ +#define icxLimitD_void ((double (*)(void *, double *))icxLimitD) /* Cast with void 1st arg */ +static double icxLimit(icxLuLut *p, double *in); /* For input */ +static int icxLuLut_init_clut_camclip(icxLuLut *p); + +/* Debug overall lookup */ +#ifdef DEBUG_OLUT +#undef DBOL +#ifdef CNDTRACE +#define DBOL(xxx) if (p->trace) printf xxx ; +#else +#define DBOL(xxx) printf xxx ; +#endif +#else +#undef DBOL +#define DBOL(xxx) +#endif + +/* Debug reverse lookup */ +#ifdef DEBUG_RLUT +#undef DBR +#ifdef CNDTRACE +#define DBR(xxx) if (p->trace) printf xxx ; +#else +#define DBR(xxx) printf xxx ; +#endif +#else +#undef DBR +#define DBR(xxx) +#endif + +/* 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 + +/* ========================================================== */ +/* xicc lookup code. */ +/* ========================================================== */ + +/* Forward and Backward Multi-Dimensional Interpolation type conversion */ +/* Return 0 on success, 1 if clipping occured, 2 on other error */ + +/* Components of overall lookup, in order */ + +int icxLuLut_in_abs(icxLuLut *p, double *out, double *in) { + int rv = 0; + + if (p->ins == icxSigJabData) { + DBOL(("xlut in_abs: CAM in = %s\n", icmPdv(p->inputChan, in))); + p->cam->cam_to_XYZ(p->cam, out, in); + DBOL(("xlut in_abs: XYZ = %s\n", icmPdv(p->inputChan, out))); + /* Hack to prevent CAM02 weirdness being amplified by inv_abs() */ + /* or 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; + DBOL(("xlut in_abs: after clipping -Y %s\n",icmPdv(p->outputChan, out))); + } + rv |= ((icmLuLut *)p->plu)->in_abs((icmLuLut *)p->plu, out, out); + DBOL(("xlut in_abs: XYZ out = %s\n", icmPdv(p->inputChan, out))); + } else { + DBOL(("xlut in_abs: PCS in = %s\n", icmPdv(p->inputChan, in))); + rv |= ((icmLuLut *)p->plu)->in_abs((icmLuLut *)p->plu, out, in); + DBOL(("xlut in_abs: PCS out = %s\n", icmPdv(p->inputChan, out))); + } + + return rv; +} + +/* Possible matrix lookup */ +/* input->input (not distinguishing matrix altered input values) */ +int icxLuLut_matrix(icxLuLut *p, double *out, double *in) { + int rv = 0; + rv |= ((icmLuLut *)p->plu)->matrix((icmLuLut *)p->plu, out, in); + return rv; +} + +/* Do input -> input' lookup */ +int icxLuLut_input(icxLuLut *p, double *out, double *in) { +#ifdef NEVER + return ((icmLuLut *)p->plu)->input((icmLuLut *)p->plu, out, in); +#else + int rv = 0; + co tc; + int i; + for (i = 0; i < p->inputChan; i++) { + tc.p[0] = in[i]; + rv |= p->inputTable[i]->interp(p->inputTable[i], &tc); + out[i] = tc.v[0]; + } + return rv; +#endif +} + +/* Do input'->output' lookup, with aux' return */ +/* (The aux' is just extracted from the in' values) */ +int icxLuLut_clut_aux(icxLuLut *p, +double *out, /* output' value */ +double *oink, /* If not NULL, return amount input is over the ink limit, 0 if not */ +double *auxv, /* If not NULL, return aux value used (packed) */ +double *in /* input' value */ +) { + int rv = 0; + co tc; + int i; + + for (i = 0; i < p->inputChan; i++) + tc.p[i] = in[i]; + rv |= p->clutTable->interp(p->clutTable, &tc); + for (i = 0; i < p->outputChan; i++) + out[i] = tc.v[i]; + + if (auxv != NULL) { + int ee = 0; + for (i = 0; i < p->clutTable->di; i++) { + double v = in[i]; + if (p->auxm[i] != 0) { + auxv[ee] = v; + ee++; + } + } + } + + if (oink != NULL) { + double lim = 0.0; + + if (p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0) { + lim = icxLimitD(p, in); + if (lim < 0.0) + lim = 0.0; + } + *oink = lim; + } + + return rv; +} + +/* Do input'->output' lookup */ +int icxLuLut_clut(icxLuLut *p, double *out, double *in) { +#ifdef NEVER + return ((icmLuLut *)p->plu)->clut((icmLuLut *)p->plu, out, in); +#else + return icxLuLut_clut_aux(p, out, NULL, NULL, in); +#endif +} + +/* Do output'->output lookup */ +int icxLuLut_output(icxLuLut *p, double *out, double *in) { + int rv = 0; + + if (p->mergeclut == 0) { +#ifdef NEVER + rv = ((icmLuLut *)p->plu)->output((icmLuLut *)p->plu, out, in); +#else + co tc; + int i; + for (i = 0; i < p->outputChan; i++) { + tc.p[0] = in[i]; + rv |= p->outputTable[i]->interp(p->outputTable[i], &tc); + out[i] = tc.v[0]; + } +#endif + } else { + int i; + for (i = 0; i < p->outputChan; i++) + out[i] = in[i]; + } + return rv; +} + +/* Relative to absolute conversion + PCS to PCS override (Effective PCS) conversion */ +int icxLuLut_out_abs(icxLuLut *p, double *out, double *in) { + int rv = 0; + if (p->mergeclut == 0) { + DBOL(("xlut out_abs: PCS in = %s\n", icmPdv(p->outputChan, in))); + + rv |= ((icmLuLut *)p->plu)->out_abs((icmLuLut *)p->plu, out, in); + + DBOL(("xlut out_abs: ABS PCS out = %s\n", icmPdv(p->outputChan, out))); + + if (p->outs == icxSigJabData) { + p->cam->XYZ_to_cam(p->cam, out, out); + + DBOL(("xlut out_abs: CAM out = %s\n", icmPdv(p->outputChan, out))); + } + } else { + int i; + for (i = 0; i < p->outputChan; i++) + out[i] = in[i]; + } + + return rv; +} + +/* Overall lookup */ +static int +icxLuLut_lookup ( +icxLuBase *pp, /* This */ +double *out, /* Vector of output values */ +double *in /* Vector of input values */ +) { + icxLuLut *p = (icxLuLut *)pp; + int rv = 0; + double temp[MAX_CHAN]; + + DBOL(("xicclu: in = %s\n", icmPdv(p->inputChan, in))); + rv |= p->in_abs (p, temp, in); + DBOL(("xicclu: after abs = %s\n", icmPdv(p->inputChan, temp))); + rv |= p->matrix (p, temp, temp); + DBOL(("xicclu: after matrix = %s\n", icmPdv(p->inputChan, temp))); + rv |= p->input (p, temp, temp); + DBOL(("xicclu: after inout = %s\n", icmPdv(p->inputChan, temp))); + rv |= p->clut (p, out, temp); + DBOL(("xicclu: after clut = %s\n", icmPdv(p->outputChan, out))); + if (p->mergeclut == 0) { + rv |= p->output (p, out, out); + DBOL(("xicclu: after output = %s\n", icmPdv(p->outputChan, out))); + rv |= p->out_abs (p, out, out); + DBOL(("xicclu: after outabs = %s\n", icmPdv(p->outputChan, 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 icxLuLut_fwd_relpcs_outpcs( +icxLuBase *pp, +icColorSpaceSignature is, /* Input space, XYZ or Lab */ +double *out, double *in) { + icxLuLut *p = (icxLuLut *)pp; + + + /* Convert to the ICC PCS */ + if (is == icSigLabData && p->natpcs == icSigXYZData) { + DBS(("fwd_relpcs_outpcs: Lab in = %s\n", icmPdv(p->inputChan, in))); + icmLab2XYZ(&icmD50, out, in); + DBS(("fwd_relpcs_outpcs: XYZ = %s\n", icmPdv(p->inputChan, out))); + } else if (is == icSigXYZData && p->natpcs == icSigLabData) { + DBS(("fwd_relpcs_outpcs: XYZ in = %s\n", icmPdv(p->inputChan, in))); + icmXYZ2Lab(&icmD50, out, in); + DBS(("fwd_relpcs_outpcs: Lab = %s\n", icmPdv(p->inputChan, out))); + } else { + DBS(("fwd_relpcs_outpcs: PCS in = %s\n", icmPdv(p->inputChan, in))); + icmCpy3(out, in); + } + + /* Convert to absolute */ + ((icmLuLut *)p->plu)->out_abs((icmLuLut *)p->plu, out, out); + + DBS(("fwd_relpcs_outpcs: abs PCS = %s\n", icmPdv(p->inputChan, out))); + + if (p->outs == icxSigJabData) { + + /* Convert to CAM */ + p->cam->XYZ_to_cam(p->cam, out, out); + + DBS(("fwd_relpcs_outpcs: Jab = %s\n", icmPdv(p->inputChan, out))); + } +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Components of inverse lookup, in order */ + +/* Utility function - compute the clip vector direction. */ +/* return NULL if vector clip isn't used. */ +double *icxClipVector( +icxClip *p, /* Clipping setup information */ +double *in, /* Target point */ +double *cdirv /* Space for returned clip vector */ +) { + int f; + if (p->nearclip != 0) + return NULL; /* Doing nearest clipping, not vector */ + + /* Default is simple vector clip */ + for (f = 0; f < p->fdi; f++) + cdirv[f] = p->ocent[f] - in[f]; /* Clip towards output gamut center */ + + if (p->ocentl != 0.0) { /* Graduated vector clip */ + double cvl, nll; + + /* Compute (negative) clip vector length */ + for (cvl = 0.0, f = 0; f < p->fdi; f++) { + cvl += cdirv[f] * cdirv[f]; + } + cvl = sqrt(cvl); + if (cvl > 1e-8) { + /* Dot product of clip vector and clip center line */ + for (nll = 0.0, f = 0; f < p->fdi; f++) + nll -= cdirv[f] * p->ocentv[f]; /* (Fix -ve clip vector) */ + nll /= (p->ocentl * p->ocentl); /* Normalised location along line */ + + /* Limit to line */ + if (nll < 0.0) + nll = 0.0; + else if (nll > 1.0) + nll = 1.0; + + if (p->LabLike) { + /* Aim more towards center for saturated targets */ + double sat = sqrt(in[1] * in[1] + in[2] * in[2]); + nll += sat/150.0 * (0.5 - nll); + } + + /* Compute target clip direction */ + for (f = 0; f < p->fdi; f++) + cdirv[f] = p->ocent[f] + nll * p->ocentv[f] - in[f]; + } + } + + return cdirv; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Typically inv is used to invert a device->PCS table, */ +/* so it is doing a PCS->device conversion. */ +/* This doesn't always have to be the case though. */ + +/* PCS override (Effective PCS) to PCS conversion + absolute to relative conversion */ +int icxLuLut_inv_out_abs(icxLuLut *p, double *out, double *in) { + int rv = 0; + + DBR(("\nicxLuLut_inv_out_abs got PCS %s\n",icmPdv(p->outputChan, in))); + + if (p->mergeclut == 0) { + if (p->outs == icxSigJabData) { + p->cam->cam_to_XYZ(p->cam, out, in); + DBR(("icxLuLut_inv_out_abs after cam2XYZ %s\n",icmPdv(p->outputChan, out))); + /* Hack to prevent CAM02 weirdness being amplified by inv_out_abs() */ + /* or 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; + DBR(("icxLuLut_inv_out_abs after clipping -Y %s\n",icmPdv(p->outputChan, out))); + } + rv |= ((icmLuLut *)p->plu)->inv_out_abs((icmLuLut *)p->plu, out, out); + DBR(("icxLuLut_inv_out_abs after icmLut inv_out_abs %s\n",icmPdv(p->outputChan, out))); + } else { + rv |= ((icmLuLut *)p->plu)->inv_out_abs((icmLuLut *)p->plu, out, in); + DBR(("icxLuLut_inv_out_abs after icmLut inv_out_abs %s\n",icmPdv(p->outputChan, out))); + } + } else { + int i; + for (i = 0; i < p->outputChan; i++) + out[i] = in[i]; + } + DBR(("icxLuLut_inv_out_abs returning PCS %f %f %f\n",out[0],out[1],out[2])) + return rv; +} + +/* Do output->output' inverse lookup */ +int icxLuLut_inv_output(icxLuLut *p, double *out, double *in) { + int rv = 0; + DBR(("icxLuLut_inv_output got PCS %f %f %f\n",in[0],in[1],in[2])) + if (p->mergeclut == 0) { +#ifdef NEVER + rv = ((icmLuLut *)p->plu)->inv_output((icmLuLut *)p->plu, out, in); +#else + int i,j; + int nsoln; /* Number of solutions found */ + co pp[MAX_INVSOLN]; /* Room for all the solutions found */ + double cdir; + + for (i = 0; i < p->outputChan; i++) { + pp[0].p[0] = p->outputClipc[i]; + pp[0].v[0] = in[i]; + cdir = p->outputClipc[i] - in[i]; /* Clip towards output range */ + + nsoln = p->outputTable[i]->rev_interp ( + p->outputTable[i], /* this */ + RSPL_NEARCLIP, /* Clip to nearest (faster than vector) */ + MAX_INVSOLN, /* Maximum number of solutions allowed for */ + NULL, /* No auxiliary input targets */ + &cdir, /* Clip vector direction and length */ + pp); /* Input and output values */ + + if (nsoln & RSPL_DIDCLIP) + rv = 1; + + nsoln &= RSPL_NOSOLNS; /* Get number of solutions */ + + if (nsoln == 1) { /* Exactly one solution */ + j = 0; + } else if (nsoln == 0) { /* Zero solutions. This is unexpected. */ + error("xlut: Unexpected failure to find reverse solution for output table"); + return 2; + } else { /* Multiple solutions */ + /* Use a simple minded resolution - choose the one closest to the center */ + double bdist = 1e300; + int bsoln = 0; + /* Don't expect this - 1D luts are meant to be monotonic */ + warning("1D lut inversion got %d reverse solutions\n",nsoln); + warning("solution 0 = %f\n",pp[0].p[0]); + warning("solution 1 = %f\n",pp[1].p[0]); + for (j = 0; j < nsoln; j++) { + double tt; + tt = pp[i].p[0] - p->outputClipc[i]; + tt *= tt; + if (tt < bdist) { /* Better solution */ + bdist = tt; + bsoln = j; + } + } + j = bsoln; + } + out[i] = pp[j].p[0]; + } + +#endif /* NEVER */ + } else { + int i; + for (i = 0; i < p->outputChan; i++) + out[i] = in[i]; + } + DBR(("icxLuLut_inv_output returning PCS' %f %f %f\n",out[0],out[1],out[2])) + return rv; +} + +/* Ink limit+gamut limit calculation function for xLuLut. */ +/* Returns < 0.0 if input value is within limits, */ +/* > 0.0 if outside limits. Value is proportinal to distance to limits. */ +/* We implement device gamut check to improve utility outside rspl, */ +/* in optimisation routines. */ +/* The limits are assumed to be post calibrated device values (ie. real */ +/* final device values) */ +static double icxLimit( +icxLuLut *p, +double *in +) { + double cin[MAX_CHAN]; /* Calibrated input values */ + double tlim, klim; + double ovr, val; + int e; + + if (p->pp->cal != NULL) { /* We have device calibration information */ + p->pp->cal->interp(p->pp->cal, cin, in); + } else { + for (e = 0; e < p->inputChan; e++) + cin[e] = in[e]; + } + + if ((tlim = p->ink.tlimit) < 0.0) + tlim = (double)p->inputChan; /* Default is no limit */ + + if ((klim = p->ink.klimit) < 0.0) + klim = 1.0; + + /* Compute amount outside total limit */ + { /* No calibration */ + double sum; + for (sum = 0.0, e = 0; e < p->inputChan; e++) + sum += cin[e]; + val = sum - tlim; + } + + /* Compute amount outside black limit */ + if (p->ink.klimit >= 0.0) { + double kval = 0.0; + switch(p->natis) { + case icSigCmykData: + kval = cin[3] - klim; + break; + default: + /* NOTE !!! p->kch isn't being initialized !!! */ + if (p->kch >= 0) { + kval = cin[p->kch] - klim; + } else { + error("xlut: Unknown colorspace when black limit specified"); + } + } + if (kval > val) + val = kval; + } + /* Compute amount outside device value limits 0.0 - 1.0 */ + for (ovr = -1.0, e = 0; e < p->inputChan; e++) { + if (in[e] < 0.0) { /* ! Not cin[] */ + if (-in[e] > ovr) + ovr = -in[e]; + } else if (in[e] > 1.0) { + if ((in[e] - 1.0) > ovr) + ovr = in[e] - 1.0; + } + } + if (ovr > val) + val = ovr; + + return val; +} + +/* Same as above, but works with input' values */ +/* (If an ink limit is being used we assume that the */ +/* input space is not PCS, hence inv_in_abs() is doing nothing) */ +static double icxLimitD( +icxLuLut *p, +double *ind +) { + double in[MAX_CHAN]; + co tc; + int e; + + /* Convert input' to input through revinput Luts (for speed) */ + for (e = 0; e < p->inputChan; e++) { + tc.p[0] = ind[e]; + p->revinputTable[e]->interp(p->revinputTable[e], &tc); + in[e] = tc.v[0]; + } + + return icxLimit(p, in); +} + +/* Ink limit+gamut limit clipping function for xLuLut (CMYK). */ +/* Return nz if there was clipping */ +static int icxDoLimit( +icxLuLut *p, +double *out, +double *in +) { + double tlim, klim = -1.0; + double sum; + int clip = 0, e; + int kch = -1; + + for (e = 0; e < p->inputChan; e++) + out[e] = in[e]; + + if ((tlim = p->ink.tlimit) < 0.0) + tlim = (double)p->inputChan; + + if ((klim = p->ink.klimit) < 0.0) + klim = 1.0; + + /* Clip black */ + if (p->natis == icSigCmykData) + kch = 3; + else + kch = p->kch; + + if (kch >= 0) { + if (out[p->kch] > klim) { + out[p->kch] = klim; + clip = 1; + } + } + + /* Compute amount outside total limit */ + for (sum = 0.0, e = 0; e < p->inputChan; e++) + sum += out[e]; + + if (sum > tlim) { + clip = 1; + sum /= (double)p->inputChan; + for (e = 0; e < p->inputChan; e++) + out[e] -= sum; + } + return clip; +} + +#ifdef NEVER +#undef DBK +#define DBK(xxx) printf xxx ; +#else +#undef DBK +#define DBK(xxx) +#endif + +/* helper function that creates our standard K locus curve value, */ +/* given the curve parameters, and the normalised L 0.0 - 1.0 value. */ +/* No filtering version. */ +/* !!! Should add K limit in here so that smoothing takes it into account !!! */ +static double icxKcurveNF(double L, icxInkCurve *c) { + double Kstpo, Kenpo, Kstle, Kenle; + double rv; + + DBK(("icxKcurve got L = %f, smth %f skew %f, Parms %f %f %f %f %f\n",L, c->Ksmth, c->Kskew, c->Kstle, c->Kstpo, c->Kenpo, c->Kenle, c->Kshap)); + + /* Invert sense of L, so that 0.0 = white, 1.0 = black */ + L = 1.0 - L; + + /* Clip L, just in case */ + if (L < 0.0) { + L = 0.0; + } else if (L > 1.0) { + L = 1.0; + } + DBK(("Clipped inverted L = %f\n",L)); + + /* Make sure breakpoints are ordered */ + if (c->Kstpo < c->Kenpo) { + Kstle = c->Kstle; + Kstpo = c->Kstpo; + Kenpo = c->Kenpo; + Kenle = c->Kenle; + } else { /* They're swapped */ + Kstle = c->Kenle; + Kstpo = c->Kenpo; + Kenpo = c->Kstpo; + Kenle = c->Kstle; + } + + if (L <= Kstpo) { + /* We are at white level */ + rv = Kstle; + DBK(("At white level %f\n",rv)); + } else if (L >= Kenpo) { + /* We are at black level */ + rv = Kenle; + DBK(("At black level %f\n",rv)); + } else { + double Lp, g; + /* We must be on the curve from start to end levels */ + + Lp = (L - Kstpo)/(Kenpo - Kstpo); + + DBK(("Curve position %f\n",Lp)); + + Lp = pow(Lp, c->Kskew); + + DBK(("Skewed curve position %f\n",Lp)); + + g = c->Kshap/2.0; + DBK(("Curve bf %f, g %g\n",Lp,g)); + + /* A value of 0.5 will be tranlated to g */ + Lp = Lp/((1.0/g - 2.0) * (1.0 - Lp) + 1.0); + + DBK(("Skewed shaped %f\n",Lp)); + + Lp = pow(Lp, 1.0/c->Kskew); + + DBK(("Shaped %f\n",Lp)); + + /* Transition between start end end levels */ + rv = Lp * (Kenle - Kstle) + Kstle; + + DBK(("Scaled to start and end levele %f\n",rv)); + } + + DBK(("Returning %f\n",rv)); + return rv; +} + +#ifdef DBK +#undef DBK +#define DBK(xxx) +#endif + + +#ifdef NEVER +#undef DBK +#define DBK(xxx) printf xxx ; +#else +#undef DBK +#define DBK(xxx) +#endif + +/* Same as above, but implement transition filters accross inflection points. */ +/* (The convolultion filter previously used could be */ +/* re-instituted if something was done about compressing */ +/* the filter at the boundaries so that the levels are met.) */ +static double icxKcurve(double L, icxInkCurve *c) { + +#ifdef DISABLE_KCURVE_FILTER + return icxKcurveNF(L, c); + +#else /* !DISABLE_KCURVE_FILTER */ + + double Kstpo, Kenpo, Kstle, Kenle, Ksmth; + double rv; + + DBK(("icxKcurve got L = %f, smth %f skew %f, Parms %f %f %f %f %f\n",L, c->Ksmth, c->Kskew, c->Kstle, c->Kstpo, c->Kenpo, c->Kenle, c->Kshap)); + + /* Invert sense of L, so that 0.0 = white, 1.0 = black */ + L = 1.0 - L; + + /* Clip L, just in case */ + if (L < 0.0) { + L = 0.0; + } else if (L > 1.0) { + L = 1.0; + } + DBK(("Clipped inverted L = %f\n",L)); + + /* Make sure breakpoints are ordered */ + if (c->Kstpo < c->Kenpo) { + Kstle = c->Kstle; + Kstpo = c->Kstpo; + Kenpo = c->Kenpo; + Kenle = c->Kenle; + } else { /* They're swapped */ + Kstle = c->Kenle; + Kstpo = c->Kenpo; + Kenpo = c->Kstpo; + Kenle = c->Kstle; + } + Ksmth = c->Ksmth; + + /* Curve value at point */ + rv = icxKcurveNF(1.0 - L, c); + + DBK(("Raw output at iL = %f, rv\n",L)); + + /* Create filtered value */ + { + double wbs, wbe; /* White transitioin start, end */ + double wbl, wfv; /* White blend factor, value at filter band */ + + double mt; /* Middle of the two transitions */ + + double bbs, bbe; /* Black transitioin start, end */ + double bbl, bfv; /* Black blend factor, value at filter band */ + + wbs = Kstpo - Ksmth; + wbe = Kstpo + Ksmth; + + bbs = Kenpo - 1.0 * Ksmth; + bbe = Kenpo + 1.0 * Ksmth; + + mt = 0.5 * (wbe + bbs); + + /* Make sure that the whit & black transition regions */ + /* don't go out of bounts or overlap */ + if (wbs < 0.0) { + wbe += wbs; + wbs = 0.0; + } + if (bbe > 1.0) { + bbs += (bbe - 1.0); + bbe = 1.0; + } + + if (wbe > mt) { + wbs += (wbe - mt); + wbe = mt; + } + + if (bbs < mt) { + bbe += (mt - bbs); + bbs = mt; + } + + DBK(("Transition windows %f - %f, %f - %f\n",wbs, wbe, bbw, bbe)); + if (wbs < wbe) { + wbl = (L - wbe)/(wbs - wbe); + + if (wbl > 0.0 && wbl < 1.0) { + wfv = icxKcurveNF(1.0 - wbe, c); + DBK(("iL = %f, wbl = %f, wfv = %f\n",L,Kstpo,wbl,wfv)); + + wbl = 1.0 - pow(1.0 - wbl, 2.0); + rv = wbl * Kstle + (1.0 - wbl) * wfv; + } + } + if (bbs < bbe) { + bbl = (L - bbe)/(bbs - bbe); + + if (bbl > 0.0 && bbl < 1.0) { + bfv = icxKcurveNF(1.0 - bbs, c); + DBK(("iL = %f, bbl = %f, bfv = %f\n",L,Kstpo,bbl,bfv)); + + bbl = pow(bbl, 2.0); + rv = bbl * bfv + (1.0 - bbl) * Kenle; + } + } + } + + /* To be safe */ + if (rv < 0.0) + rv = 0.0; + else if (rv > 1.0) + rv = 1.0; + + DBK(("Returning %f\n",rv)); + return rv; +#endif /* !DISABLE_KCURVE_FILTER */ +} + +#ifdef DBK +#undef DBK +#define DBK(xxx) +#endif + +/* Do output'->input' lookup with aux details. */ +/* Note that out[] will be used as the inking value if icxKrule is */ +/* icxKvalue, icxKlocus, icxKl5l or icxKl5lk, and that the auxiliar values, PCS ranges */ +/* and icxKrule value will all be evaluated in output->input space (not ' space). */ +/* Note that the ink limit will be computed after converting input' to input, auxt */ +/* will override the inking rule, and auxr[] reflects the available auxiliary range */ +/* that the locus was to choose from, and auxv[] was the actual auxiliary used. */ +/* Returns clip status. */ +int icxLuLut_inv_clut_aux( +icxLuLut *p, +double *out, /* Function return values, plus aux value or locus target input if auxt == NULL */ +double *auxv, /* If not NULL, return aux value used (packed) */ +double *auxr, /* If not NULL, return aux locus range (packed, 2 at a time) */ +double *auxt, /* If not NULL, specify the aux target for this lookup (override ink) */ +double *clipd, /* If not NULL, return DE to gamut on clipi, 0 for not clip */ +double *in /* Function input values to invert (== clut output' values) */ +) { + co pp[MAX_INVSOLN]; /* Room for all the solutions found */ + int nsoln; /* Number of solutions found */ + double *cdir, cdirv[MXDO]; /* Clip vector direction and length */ + int e,f,i; + int fdi = p->clutTable->fdi; + int flags = 0; /* reverse interp flags */ + int xflags = 0; /* extra clip/noclip flags */ + double tin[MXDO]; /* PCS value to be inverted */ + double cdist = 0.0; /* clip DE */ + int crv = 0; /* Return value - set to 1 if clipped */ + + if (p->nearclip != 0) + flags |= RSPL_NEARCLIP; /* Use nearest clipping rather than clip vector */ + + DBR(("inv_clut_aux input is %f %f %f\n",in[0], in[1], in[2])) + + if (auxr != NULL) { /* Set a default locus range */ + int ee = 0; + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + auxr[ee++] = 1e60; + auxr[ee++] = -1e60; + } + } + } + + /* Setup for reverse lookup */ + for (f = 0; f < fdi; f++) + pp[0].v[f] = in[f]; /* Target value */ + + /* Compute clip vector, if any */ + cdir = icxClipVector(&p->clip, in, cdirv); + + if (p->clutTable->di > fdi) { /* ie. CMYK->Lab, there will be ambiguity */ + double min[MXDI], max[MXDI]; /* Auxiliary locus range */ + +#ifdef REPORT_LOCUS_SEGMENTS /* Examine how many segments there are */ + { /* ie. CMYK->Lab, there will be ambiguity */ + double smin[10][MXRI], smax[10][MXRI]; /* Auxiliary locus segment ranges */ + double min[MXRI], max[MXRI]; /* Auxiliary min and max locus range */ + + nsoln = p->clutTable->rev_locus_segs( + p->clutTable, /* rspl object */ + p->auxm, /* Auxiliary mask */ + pp, /* Input target and output solutions */ + 10, /* Maximum number of solutions to return */ + smin, smax); /* Returned locus of valid auxiliary values */ + + if (nsoln != 0) { + /* Convert the locuses from input' -> input space */ + /* and get overall min/max locus range */ + for (e = 0; e < p->clutTable->di; e++) { + co tc; + /* (Is speed more important than precision ?) */ + if (p->auxm[e] != 0) { + for (i = 0; i < nsoln; i++) { + tc.p[0] = smin[i][e]; + p->revinputTable[e]->interp(p->revinputTable[e], &tc); + smin[i][e] = tc.v[0]; + tc.p[0] = smax[i][e]; + p->revinputTable[e]->interp(p->revinputTable[e], &tc); + smax[i][e] = tc.v[0]; + printf(" Locus seg %d:[%d] %f -> %f\n",i, e, smin[i][e], smax[i][e]); + } + } + } + } + } +#endif /* REPORT_LOCUS_SEGMENTS */ + + /* Compute auxiliary locus on the fly. This is in dev' == input' space. */ + nsoln = p->clutTable->rev_locus( + p->clutTable, /* rspl object */ + p->auxm, /* Auxiliary mask */ + pp, /* Input target and output solutions */ + min, max); /* Returned locus of valid auxiliary values */ + + if (nsoln == 0) { + xflags |= RSPL_WILLCLIP; /* No valid locus, so we expect to have to clip */ +#ifdef DEBUG_RLUT + printf("inv_clut_aux: no valid locus, expect clip\n"); +#endif + /* Make sure that the auxiliar value is initialized, */ + /* even though it won't have any effect. */ + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + pp[0].p[e] = 0.5; + } + } + + } else { /* Got a valid locus */ + + /* Convert the locuses from input' -> input space */ + for (e = 0; e < p->clutTable->di; e++) { + co tc; + /* (Is speed more important than precision ?) */ + if (p->auxm[e] != 0) { + tc.p[0] = min[e]; + p->revinputTable[e]->interp(p->revinputTable[e], &tc); + min[e] = tc.v[0]; + tc.p[0] = max[e]; + p->revinputTable[e]->interp(p->revinputTable[e], &tc); + max[e] = tc.v[0]; + } + } + + if (auxr != NULL) { /* Report the locus range */ + int ee = 0; + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + auxr[ee++] = min[e]; + auxr[ee++] = max[e]; + } + } + } + + if (auxt != NULL) { /* overiding auxiliary target */ + int ee = 0; + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + double iv = auxt[ee++]; + if (iv < min[e]) + iv = min[e]; + else if (iv > max[e]) + iv = max[e]; + pp[0].p[e] = iv; + } + } + DBR(("inv_clut_aux: aux %f from auxt[] %f\n",pp[0].p[3],auxt[0])) + } else if (p->ink.k_rule == icxKvalue) { + /* Implement the auxiliary inking rule */ + /* Target auxiliary values are provided in out[] K value */ + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + double iv = out[e]; /* out[] holds aux target value */ + if (iv < min[e]) + iv = min[e]; + else if (iv > max[e]) + iv = max[e]; + pp[0].p[e] = iv; + } + } + DBR(("inv_clut_aux: aux %f from out[0] K target %f min %f max %f\n",pp[0].p[3],out[3],min[3],max[3])) + } else if (p->ink.k_rule == icxKlocus) { + /* Set target auxliary input values from values in out[] and locus */ + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + double ii, iv; + ii = out[e]; /* Input ink locus */ + iv = min[e] + ii * (max[e] - min[e]); /* Output ink from locus */ + if (iv < min[e]) + iv = min[e]; + else if (iv > max[e]) + iv = max[e]; + pp[0].p[e] = iv; + } + } + DBR(("inv_clut_aux: aux %f from out[0] locus %f min %f max %f\n",pp[0].p[3],out[3],min[3],max[3])) + } else { /* p->ink.k_rule == icxKluma5 || icxKluma5k || icxKl5l || icxKl5lk */ + /* Auxiliaries are driven by a rule and the output values */ + double rv, L; + + /* If we've got a mergeclut, then the PCS' is the same as the */ + /* effective PCS, and we need to convert to native PCS */ + if (p->mergeclut) { + p->mergeclut = 0; /* Hack to be able to use inv_out_abs() */ + icxLuLut_inv_out_abs(p, tin, in); + p->mergeclut = 1; + + } else { + /* Convert native PCS' to native PCS values */ + p->output(p, tin, in); + } + + /* Figure out Luminance number */ + if (p->natos == icSigXYZData) { + icmXYZ2Lab(&icmD50, tin, tin); + } else if (p->natos != icSigLabData) { /* Hmm. that's unexpected */ + error("Assert: xlut K locus, unexpected native pcs of 0x%x\n",p->natos); + } + L = 0.01 * tin[0]; + DBR(("inv_clut_aux: aux from Luminance, raw L = %f\n",L)); + + /* Normalise L to its possible range from min to max */ + L = (L - p->Lmin)/(p->Lmax - p->Lmin); + DBR(("inv_clut_aux: Normalize L = %f\n",L)); + + /* Convert L to curve value */ + rv = icxKcurve(L, &p->ink.c); + DBR(("inv_clut_aux: icxKurve lookup returns = %f\n",rv)); + + if (p->ink.k_rule == icxKluma5) { /* Curve is locus value */ + + /* Set target black as K fraction within locus */ + + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + pp[0].p[e] = min[e] + rv * (max[e] - min[e]); + } + } + DBR(("inv_clut_aux: aux %f from locus %f min %f max %f\n",pp[0].p[3],rv,min[3],max[3])) + + } else if (p->ink.k_rule == icxKluma5k) { /* Curve is K value */ + + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + double iv = rv; + if (iv < min[e]) /* Clip to locus */ + iv = min[e]; + else if (iv > max[e]) + iv = max[e]; + pp[0].p[e] = iv; + } + } + DBR(("inv_clut_aux: aux %f from out[0] K target %f min %f max %f\n",pp[0].p[3],rv,min[3],max[3])) + + } else { /* icxKl5l || icxKl5lk */ + /* Create second curve, and use input locus to */ + /* blend between */ + + double rv2; /* Upper limit */ + + /* Convert L to max curve value */ + rv2 = icxKcurve(L, &p->ink.x); + + if (rv2 < rv) { /* Ooops - better swap. */ + double tt; + tt = rv; + rv = rv2; + rv2 = tt; + } + + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + if (p->ink.k_rule == icxKl5l) { + double ii; + ii = out[e]; /* Input K locus */ + if (ii < 0.0) + ii = 0.0; + else if (ii > 1.0) + ii = 1.0; + ii = (1.0 - ii) * rv + ii * rv2;/* Blend between locus rule curves */ + /* Out ink from output locus */ + pp[0].p[e] = min[e] + ii * (max[e] - min[e]); + } else { + double iv; + iv = out[e]; /* Input K level */ + if (iv < rv) /* Constrain to curves */ + iv = rv; + else if (iv > rv2) + iv = rv2; + pp[0].p[e] = iv; + } + } + } + DBR(("inv_clut_aux: aux %f from 2 curves\n",pp[0].p[3])) + } + } + + /* Convert to input/dev aux target to input'/dev' space for rspl inversion */ + for (e = 0; e < p->clutTable->di; e++) { + double tv, bv = 0.0, bd = 1e6; + co tc; + if (p->auxm[e] != 0) { + tv = pp[0].p[e]; + /* Clip to overall locus range (belt and braces) */ + if (tv < min[e]) + tv = min[e]; + if (tv > max[e]) + tv = max[e]; + tc.p[0] = tv; + p->inputTable[e]->interp(p->inputTable[e], &tc); + pp[0].p[e] = tc.v[0]; + } + } + + xflags |= RSPL_EXACTAUX; /* Since we confine aux to locus */ + +#ifdef DEBUG_RLUT + printf("inv_clut_aux computed aux values "); + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) + printf("%d: %f ",e,pp[0].p[e]); + } + printf("\n"); +#endif /* DEBUG_RLUT */ + } + + if (clipd != NULL) { /* Copy pp.v[] to compute clip DE */ + for (f = 0; f < fdi; f++) + tin[f] = pp[0].v[f]; + } + + /* Find reverse solution with target auxiliaries */ + /* We choose the closest aux at or above the target */ + /* to try and avoid glitches near black due to */ + /* possible forked black locuses. */ + nsoln = p->clutTable->rev_interp( + p->clutTable, /* rspl object */ + RSPL_MAXAUX | flags | xflags, /* Combine all the flags */ + MAX_INVSOLN, /* Maxumum solutions to return */ + p->auxm, /* Auxiliary input chanel mask */ + cdir, /* Clip vector direction and length */ + pp); /* Input target and output solutions */ + /* returned solutions in pp[0..retval-1].p[] */ + + } else { + DBR(("inv_clut_aux needs no aux value\n")) + + if (clipd != NULL) { /* Copy pp.v[] to compute clip DE */ + for (f = 0; f < fdi; f++) + tin[f] = pp[0].v[f]; + } + + /* Color spaces don't need auxiliaries to choose from solution locus */ + nsoln = p->clutTable->rev_interp( + p->clutTable, /* rspl object */ + flags, /* No extra flags */ + MAX_INVSOLN, /* Maxumum solutions to return */ + NULL, /* No auxiliary input targets */ + cdir, /* Clip vector direction and length */ + pp); /* Input target and output solutions */ + /* returned solutions in pp[0..retval-1].p[] */ + } + if (nsoln & RSPL_DIDCLIP) + crv = 1; /* Clipped on PCS inverse lookup */ + + if (crv && clipd != NULL) { + + /* Compute the clip DE */ + for (cdist = 0.0, f = 0; f < fdi; f++) { + double tt; + tt = pp[0].v[f] - tin[f]; + cdist += tt * tt; + } + cdist = sqrt(cdist); + } + + nsoln &= RSPL_NOSOLNS; /* Get number of solutions */ + + DBR(("inv_clut_aux got %d rev_interp solutions, clipflag = %d\n",nsoln,crv)) + + /* If we clipped and we should clip in CAM Jab space, compute reverse */ + /* clip solution using our additional CAM space. */ + /* (Note that we don't support vector clip in CAM space at the moment) */ + if (crv != 0 && p->camclip && p->nearclip) { + co cpp; /* Alternate CAM space solution */ + double bf; /* Blend factor */ + + DBR(("inv_clut_aux got clip, compute CAM clip\n")) + + if (nsoln != 1) { /* This would be unexpected */ + error("Unexpected failure to return 1 solution on clip for input to output table"); + } + + if (p->cclutTable == NULL) { /* we haven't created this yet, so do so */ + if (icxLuLut_init_clut_camclip(p)) + error("Creating CAM rspl for camclip failed"); + } + + /* Setup for reverse lookup */ + DBR(("inv_clut_aux cam clip PCS in %f %f %f\n",in[0],in[1],in[2])) + + /* Convert from PCS' to (XYZ) PCS */ + ((icmLuLut *)p->absxyzlu)->output((icmLuLut *)p->absxyzlu, tin, in); + DBR(("inv_clut_aux cam clip PCS' -> PCS %f %f %f\n",tin[0],tin[1],tin[2])) + + ((icmLuLut *)p->absxyzlu)->out_abs((icmLuLut *)p->absxyzlu, tin, tin); + DBR(("inv_clut_aux cam clip abs XYZ PCS %f %f %f\n",tin[0],tin[1],tin[2])) + + p->cam->XYZ_to_cam(p->cam, tin, tin); + DBR(("inv_clut_aux cam clip PCS after XYZtoCAM %f %f %f\n",tin[0],tin[1],tin[2])) + + for (f = 0; f < fdi; f++) /* Transfer CAM targ */ + cpp.v[f] = tin[f]; + + /* Make sure that the auxiliar value is initialized, */ + /* even though it shouldn't have any effect, since should clipp. */ + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + cpp.p[e] = 0.5; + } + } + + if (p->clutTable->di > fdi) { /* ie. CMYK->Lab, there will be ambiguity */ + + nsoln = p->cclutTable->rev_interp( + p->cclutTable, /* rspl object */ + flags | xflags | RSPL_WILLCLIP, /* Combine all the flags + clip ?? */ + 1, /* Maximum solutions to return */ + p->auxm, /* Auxiliary input chanel mask */ + cdir, /* Clip vector direction and length */ + &cpp); /* Input target and output solutions */ + + } else { + nsoln = p->cclutTable->rev_interp( + p->cclutTable, /* rspl object */ + flags | RSPL_WILLCLIP, /* Because we know it will clip ?? */ + 1, /* Maximum solutions to return */ + NULL, /* No auxiliary input targets */ + cdir, /* Clip vector direction and length */ + &cpp); /* Input target and output solutions */ + } + + nsoln &= RSPL_NOSOLNS; /* Get number of solutions */ + + if (nsoln != 1) { /* This would be unexpected */ + error("Unexpected failure to return 1 solution on CAM clip for input to output table"); + } + + /* Compute the CAM clip distances */ + for (cdist = 0.0, f = 0; f < fdi; f++) { + double tt; + tt = cpp.v[f] - tin[f]; + cdist += tt * tt; + } + cdist = sqrt(cdist); + + /* Use magic number to set blend distance, and compute a blend factor. */ + /* Blend over 1 delta E, otherwise the Lab & CAM02 divergence can result */ + /* in reversals. */ + bf = cdist/CAMCLIPTRANS; /* 0.0 for PCS result, 1.0 for CAM result */ + if (bf > 1.0) + bf = 1.0; +#ifdef USECAMCLIPSPLINE + bf = bf * bf * (3.0 - 2.0 * bf); /* Convert to spline blend */ +#endif + DBR(("cdist %f, spline blend %f\n",cdist,bf)) + + /* Blend between solution values for PCS and CAM clip result. */ + /* We're hoping that the solutions are close, and expect them to be */ + /* that way when we're close to the gamut surface (since the cell */ + /* vertex values should be exact, irrespective of which space they're */ + /* in), but weird things could happen away from the surface. Weird */ + /* things can happen anyway with "clip to nearest", since this is not */ + /* guaranteed to be a smooth function, depending on the gamut surface */ + /* geometry, without taking some precaution such as clipping to a */ + /* convex hull "wrapper" to create a clip vector, which we're not */ + /* current doing. */ + DBR(("Clip blend between device:\n")) + DBR(("Lab: %f %f %f and\n",pp[0].p[0], pp[0].p[1], pp[0].p[2])) + DBR(("Jab: %f %f %f\n",cpp.p[0], cpp.p[1], cpp.p[2])) + + for (e = 0; e < p->clutTable->di; e++) { + out[e] = (1.0 - bf) * pp[0].p[e] + bf * cpp.p[e]; + } + + /* Not CAM clip case */ + } else { + + if (nsoln == 1) { /* Exactly one solution */ + i = 0; + } else if (nsoln == 0) { /* Zero solutions. This is unexpected. */ + double in_v[MXDO]; + p->output(p, in_v, pp[0].v); /* Get ICC inverse input values */ + p->out_abs(p, in_v, in_v); + error("Unexpected failure to find reverse solution for input to output table for value %f %f %f (ICC input %f %f %f)",pp[0].v[0],pp[0].v[1],pp[0].v[2], in_v[0], in_v[1], in_v[2]); + return 2; + } else { /* Multiple solutions */ + /* Use a simple minded resolution - choose the one closest to the center */ + double bdist = 1e300; + int bsoln = 0; + DBR(("got multiple reverse solutions\n")); + for (i = 0; i < nsoln; i++) { + double ss; + + for (ss = 0.0, e = 0; e < p->clutTable->di; e++) { + double tt; + tt = pp[i].p[e] - p->licent[e]; + tt *= tt; + if (tt < bdist) { /* Better solution */ + bdist = tt; + bsoln = i; + } + } + } +//printf("~1 chose %d\n",bsoln); + i = bsoln; + } + for (e = 0; e < p->clutTable->di; e++) { + /* Save solution as atractor for next one, on the basis */ + /* that it might have better continuity given pesudo-hilbert inversion path. */ + p->licent[e] = out[e] = pp[i].p[e]; /* Solution */ + } + } + + /* Sanitise auxiliary locus range and auxiliary value return */ + if (auxr != NULL || auxv != NULL) { + int ee = 0; + for (e = 0; e < p->clutTable->di; e++) { + double v = out[e]; /* Solution */ + if (p->auxm[e] != 0) { + if (auxr != NULL) { /* Make sure locus encloses actual value */ + if (auxr[2 * ee] > v) + auxr[2 * ee] = v; + if (auxr[2 * ee + 1] < v) + auxr[2 * ee + 1] = v; + } + if (auxv != NULL) { + auxv[ee] = v; + } + ee++; + } + } + } + +#ifdef CHECK_ILIMIT /* Do sanity checks on meeting ink limit */ +if (p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0) { + double sum = icxLimitD(p, out); + if (sum > 0.0) + printf("xlut assert%s: icxLuLut_inv_clut returned outside limits by %f > tlimit %f\n",crv ? " (clip)" : "", sum, p->ink.tlimit); +} +#endif + + if (clipd != NULL) { + *clipd = cdist; + DBR(("inv_clut_aux returning clip DE %f\n",cdist)) + } + + DBR(("inv_clut_aux returning %f %f %f %f\n",out[0],out[1],out[2],out[3])) + return crv; +} + +/* Do output'->input' lookup, simple version */ +/* Note than out[] will carry inking value if icxKrule is icxKvalue of icxKlocus */ +/* and that the icxKrule value will be in the input (NOT input') space. */ +/* Note that the ink limit will be computed after converting input' to input */ +int icxLuLut_inv_clut(icxLuLut *p, double *out, double *in) { + return icxLuLut_inv_clut_aux(p, out, NULL, NULL, NULL, NULL, in); +} + +/* Given the proposed auxiliary input values in in[di], */ +/* and the target output' (ie. PCS') values in out[fdi], */ +/* return the auxiliary input (NOT input' space) values as a proportion of their */ +/* possible locus in locus[di]. */ +/* This is generally used on a source CMYK profile to convey the black intent */ +/* to destination CMYK profile. */ +int icxLuLut_clut_aux_locus(icxLuLut *p, double *locus, double *out, double *in) { + co pp[1]; /* Room for all the solutions found */ + int nsoln; /* Number of solutions found */ + int e,f; + + if (p->clutTable->di > p->clutTable->fdi) { /* ie. CMYK->Lab, there will be ambiguity */ + double min[MXDI], max[MXDI]; /* Auxiliary locus range */ + + /* Setup for auxiliary locus lookup */ + for (f = 0; f < p->clutTable->fdi; f++) { + pp[0].v[f] = out[f]; /* Target output' (i.e. PCS) value */ + } + + /* Compute auxiliary locus */ + nsoln = p->clutTable->rev_locus( + p->clutTable, /* rspl object */ + p->auxm, /* Auxiliary mask */ + pp, /* Input target and output solutions */ + min, max); /* Returned locus of valid auxiliary values */ + + if (nsoln == 0) { + for (e = 0; e < p->clutTable->di; e++) + locus[e] = 0.0; /* Return some safe values */ + } else { /* Got a valid locus */ + + /* Convert the locus from input' -> input space */ + for (e = 0; e < p->clutTable->di; e++) { + co tc; + /* (Is speed more important than precision ?) */ + if (p->auxm[e] != 0) { + tc.p[0] = min[e]; + p->revinputTable[e]->interp(p->revinputTable[e], &tc); + min[e] = tc.v[0]; + tc.p[0] = max[e]; + p->revinputTable[e]->interp(p->revinputTable[e], &tc); + max[e] = tc.v[0]; + } + } + + /* Figure out the proportion of the locus */ + for (e = 0; e < p->clutTable->di; e++) { + if (p->auxm[e] != 0) { + double iv = in[e]; + if (iv <= min[e]) + locus[e] = 0.0; + else if (iv >= max[e]) + locus[e] = 1.0; + else { + double lpl = max[e] - min[e]; /* Locus path length */ + if (lpl > 1e-6) + locus[e] = (iv - min[e])/lpl; + else + locus[e] = 0.0; + } + } + } + } + } else { + /* There should be no auxiliaries */ + for (e = 0; e < p->clutTable->di; e++) + locus[e] = 0.0; /* Return some safe values */ + } + return 0; +} + +/* Do input' -> input inverse lookup */ +int icxLuLut_inv_input(icxLuLut *p, double *out, double *in) { +#ifdef NEVER + return ((icmLuLut *)p->plu)->inv_input((icmLuLut *)p->plu, out, in); +#else + int rv = 0; + int i,j; + int nsoln; /* Number of solutions found */ + co pp[MAX_INVSOLN]; /* Room for all the solutions found */ + double cdir; + + DBR(("inv_input got DEV' %f %f %f %f\n",in[0],in[1],in[2],in[3])) + + for (i = 0; i < p->inputChan; i++) { + pp[0].p[0] = p->inputClipc[i]; + pp[0].v[0] = in[i]; + cdir = p->inputClipc[i] - in[i]; /* Clip towards output range */ + + nsoln = p->inputTable[i]->rev_interp ( + p->inputTable[i], /* this */ + RSPL_NEARCLIP, /* Clip to nearest (faster than vector) */ + MAX_INVSOLN, /* Maximum number of solutions allowed for */ + NULL, /* No auxiliary input targets */ + &cdir, /* Clip vector direction and length */ + pp); /* Input and output values */ + + if (nsoln & RSPL_DIDCLIP) + rv = 1; + + nsoln &= RSPL_NOSOLNS; /* Get number of solutions */ + + if (nsoln == 1) { /* Exactly one solution */ + j = 0; + } else if (nsoln == 0) { /* Zero solutions. This is unexpected. */ + error("Unexpected failure to find reverse solution for input table"); + return 2; + } else { /* Multiple solutions */ + /* Use a simple minded resolution - choose the one closest to the center */ + double bdist = 1e300; + int bsoln = 0; + /* Don't expect this - 1D luts are meant to be monotonic */ + warning("1D lut inversion got %d reverse solutions\n",nsoln); + warning("solution 0 = %f\n",pp[0].p[0]); + warning("solution 1 = %f\n",pp[1].p[0]); + for (j = 0; j < nsoln; j++) { + double tt; + tt = pp[i].p[0] - p->inputClipc[i]; + tt *= tt; + if (tt < bdist) { /* Better solution */ + bdist = tt; + bsoln = j; + } + } + j = bsoln; + } + out[i] = pp[j].p[0]; + } + + DBR(("inv_input returning DEV %f %f %f %f\n",out[0],out[1],out[2],out[3])) + return rv; +#endif /* NEVER */ +} + +/* Possible inverse matrix lookup */ +/* (Will do nothing if input is device space) */ +int icxLuLut_inv_matrix(icxLuLut *p, double *out, double *in) { + int rv = 0; + rv |= ((icmLuLut *)p->plu)->inv_matrix((icmLuLut *)p->plu, out, in); + return rv; +} + +/* Inverse input absolute intent conversion */ +/* (Will do nothing if input is device space) */ +int icxLuLut_inv_in_abs(icxLuLut *p, double *out, double *in) { + int rv = 0; + rv |= ((icmLuLut *)p->plu)->inv_in_abs((icmLuLut *)p->plu, out, in); + + if (p->ins == icxSigJabData) { + p->cam->XYZ_to_cam(p->cam, out, out); + } + + return rv; +} + +/* Overall inverse lookup */ +/* Note that all auxiliary values are in input (NOT input') space */ +static int +icxLuLut_inv_lookup( +icxLuBase *pp, /* This */ +double *out, /* Vector of output values/input auxiliary values */ +double *in /* Vector of input values */ +) { + icxLuLut *p = (icxLuLut *)pp; + int rv = 0; + int i; + double temp[MAX_CHAN]; + + DBOL(("xiccilu: input = %s\n", icmPdv(p->outputChan, in))); + if (p->mergeclut == 0) { /* Do this if it's not merger with clut */ + rv |= p->inv_out_abs (p, temp, in); + DBOL(("xiccilu: after inv abs = %s\n", icmPdv(p->outputChan, temp))); + rv |= p->inv_output (p, temp, temp); + DBOL(("xiccilu: after inv out = %s\n", icmPdv(p->outputChan, temp))); + } else { + for (i = 0; i < p->outputChan; i++) + temp[i] = in[i]; + } + DBOL(("xiccilu: aux targ = %s\n", icmPdv(p->inputChan,out))); + rv |= p->inv_clut (p, out, temp); + DBOL(("xiccilu: after inv clut = %s\n", icmPdv(p->inputChan,out))); + rv |= p->inv_input (p, out, out); + DBOL(("xiccilu: after inv input = %s\n", icmPdv(p->inputChan,out))); + rv |= p->inv_matrix (p, out, out); + DBOL(("xiccilu: after inv matrix = %s\n", icmPdv(p->inputChan,out))); + rv |= p->inv_in_abs (p, out, out); + DBOL(("xiccilu: after inv abs = %s\n", icmPdv(p->inputChan,out))); + return rv; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* 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 icxLuLut_bwd_outpcs_relpcs( +icxLuBase *pp, +icColorSpaceSignature os, /* Output space, XYZ or Lab */ +double *out, double *in) { + icxLuLut *p = (icxLuLut *)pp; + + if (p->outs == icxSigJabData) { + DBS(("bwd_outpcs_relpcs: Jab in = %s\n", icmPdv(3, in))); + p->cam->cam_to_XYZ(p->cam, out, in); + DBS(("bwd_outpcs_relpcs: abs XYZ = %s\n", icmPdv(3, out))); + /* Hack to prevent CAM02 weirdness being amplified by */ + /* 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; + DBS(("bwd_outpcs_relpcs: after clipping -Y %s\n",icmPdv(p->outputChan, out))); + } + } else { + DBS(("bwd_outpcs_relpcs: abs PCS in = %s\n", icmPdv(3, out))); + icmCpy3(out, in); + } + + ((icmLuLut *)p->plu)->inv_out_abs((icmLuLut *)p->plu, out, out); + DBS(("bwd_outpcs_relpcs: rel PCS = %s\n", icmPdv(3, out))); + + if (os == icSigXYZData && p->natpcs == icSigLabData) { + icmLab2XYZ(&icmD50, out, out); + DBS(("bwd_outpcs_relpcs: rel XYZ = %s\n", icmPdv(3, out))); + } else if (os == icSigXYZData && p->natpcs == icSigLabData) { + icmXYZ2Lab(&icmD50, out, out); + DBS(("bwd_outpcs_relpcs: rel Lab = %s\n", icmPdv(3, out))); + } +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/* Return LuLut information */ +static void icxLuLut_get_info( + icxLuLut *p, /* this */ + icmLut **lutp, /* Pointer to icc lut type return value */ + icmXYZNumber *pcswhtp, /* Pointer to profile PCS white point return value */ + icmXYZNumber *whitep, /* Pointer to profile absolute white point return value */ + icmXYZNumber *blackp /* Pointer to profile absolute black point return value */ +) { + ((icmLuLut *)p->plu)->get_info((icmLuLut *)p->plu, lutp, pcswhtp, whitep, blackp); +} + +/* Return the underlying Lut matrix */ +static void +icxLuLut_get_matrix ( + icxLuLut *p, + double m[3][3] +) { + ((icmLuLut *)p->plu)->get_matrix((icmLuLut *)p->plu, m); +} + +static void +icxLuLut_free( +icxLuBase *pp +) { + icxLuLut *p = (icxLuLut *)pp; + int i; + + for (i = 0; i < p->inputChan; i++) { + if (p->inputTable[i] != NULL) + p->inputTable[i]->del(p->inputTable[i]); + if (p->revinputTable[i] != NULL) + p->revinputTable[i]->del(p->revinputTable[i]); + } + + if (p->clutTable != NULL) + p->clutTable->del(p->clutTable); + + if (p->cclutTable != NULL) + p->cclutTable->del(p->cclutTable); + + for (i = 0; i < p->outputChan; i++) { + if (p->outputTable[i] != NULL) + p->outputTable[i]->del(p->outputTable[i]); + } + + if (p->plu != NULL) + p->plu->del(p->plu); + + if (p->cam != NULL) + p->cam->del(p->cam); + + if (p->absxyzlu != NULL) + p->absxyzlu->del(p->absxyzlu); + + free(p); +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +static gamut *icxLuLutGamut(icxLuBase *plu, double detail); + +/* Do the basic icxLuLut creation and initialisation */ +static icxLuLut * +alloc_icxLuLut( + xicc *xicp, + icmLuBase *plu, /* Pointer to Lu we are expanding (ours) */ + int flags /* clip, merge flags */ +) { + icxLuLut *p; /* Object being created */ + icmLuLut *luluto = (icmLuLut *)plu; /* Lookup Lut type object */ + + if ((p = (icxLuLut *) calloc(1,sizeof(icxLuLut))) == NULL) + return NULL; + + p->pp = xicp; + p->plu = plu; + p->del = icxLuLut_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 = icxLuLutGamut; + p->fwd_relpcs_outpcs = icxLuLut_fwd_relpcs_outpcs; + p->bwd_outpcs_relpcs = icxLuLut_bwd_outpcs_relpcs; + p->nearclip = 0; /* Set flag defaults */ + p->mergeclut = 0; + p->noisluts = 0; + p->noipluts = 0; + p->nooluts = 0; + p->intsep = 0; + + p->lookup = icxLuLut_lookup; + p->in_abs = icxLuLut_in_abs; + p->matrix = icxLuLut_matrix; + p->input = icxLuLut_input; + p->clut = icxLuLut_clut; + p->clut_aux = icxLuLut_clut_aux; + p->output = icxLuLut_output; + p->out_abs = icxLuLut_out_abs; + + p->inv_lookup = icxLuLut_inv_lookup; + p->inv_in_abs = icxLuLut_inv_in_abs; + p->inv_matrix = icxLuLut_inv_matrix; + p->inv_input = icxLuLut_inv_input; + p->inv_clut = icxLuLut_inv_clut; + p->inv_clut_aux = icxLuLut_inv_clut_aux; + p->inv_output = icxLuLut_inv_output; + p->inv_out_abs = icxLuLut_inv_out_abs; + + p->clut_locus = icxLuLut_clut_aux_locus; + + p->get_info = icxLuLut_get_info; + p->get_matrix = icxLuLut_get_matrix; + + /* Setup all the rspl analogs of the icc Lut */ + /* NOTE: We assume that none of this relies on the flag settings, */ + /* since they will be set on our return. */ + + /* Get details of underlying, native icc color space */ + p->plu->lutspaces(p->plu, &p->natis, NULL, &p->natos, NULL, &p->natpcs); + + /* Get other details of conversion */ + p->plu->spaces(p->plu, NULL, &p->inputChan, NULL, &p->outputChan, NULL, NULL, NULL, NULL, NULL); + + /* Remember flags */ + p->flags = flags; + + /* Sanity check */ + if (p->inputChan > MXDI) { + sprintf(p->pp->err,"xicc can only handle input channels of %d or less",MXDI); + p->inputChan = MXDI; /* Avoid going outside array bounds */ + p->pp->errc = 1; + p->del((icxLuBase *)p); + return NULL; + } + if (p->outputChan > MXDO) { + sprintf(p->pp->err,"xicc can only handle output channels of %d or less",MXDO); + p->outputChan = MXDO; /* Avoid going outside array bounds */ + p->pp->errc = 1; + p->del((icxLuBase *)p); + return NULL; + } + + /* Get pointer to icmLut */ + luluto->get_info(luluto, &p->lut, NULL, NULL, NULL); + + return p; +} + +/* Initialise the clut ink limiting and black */ +/* generation information. */ +/* return 0 or error code */ +static int +setup_ink_icxLuLut( +icxLuLut *p, /* Object being initialised */ +icxInk *ink, /* inking details (NULL for default) */ +int setLminmax /* Figure the L locus for inking rule */ +) { + int devchan = p->func == icmFwd ? p->inputChan : p->outputChan; + + if (ink) { + p->ink = *ink; /* Copy the structure */ + } else { + p->ink.tlimit = 3.0; /* default ink limit of 300% */ + p->ink.klimit = -1.0; /* default no black limit */ + p->ink.KonlyLmin = 0; /* default don't use K only black as Locus Lmin */ + p->ink.k_rule = icxKluma5; /* default K generation rule */ + p->ink.c.Ksmth = ICXINKDEFSMTH; /* Default smoothing */ + p->ink.c.Kskew = ICXINKDEFSKEW; /* Default shape skew */ + p->ink.c.Kstle = 0.0; /* Min K at white end */ + p->ink.c.Kstpo = 0.0; /* Start of transition is at white */ + p->ink.c.Kenle = 1.0; /* Max K at black end */ + p->ink.c.Kenpo = 1.0; /* End transition at black */ + p->ink.c.Kshap = 1.0; /* Linear transition */ + } + + /* Normalise total and black ink limits */ + if (p->ink.tlimit <= 1e-4 || p->ink.tlimit >= (double)devchan) + p->ink.tlimit = -1.0; /* Turn ink limit off if not effective */ + if (devchan < 4 || p->ink.klimit < 0.0 || p->ink.klimit >= 1.0) + p->ink.klimit = -1.0; /* Turn black limit off if not effective */ + + /* Set the ink limit information for any reverse interpolation. */ + /* Calling this will clear the reverse interpolaton cache. */ + p->clutTable->rev_set_limit( + p->clutTable, /* this */ + p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0 ? icxLimitD_void : NULL, + /* Optional input space limit() function. */ + /* Function should evaluate in[0..di-1], and return */ + /* number that is not to exceed 0.0. NULL if not used. */ + (void *)p, /* Context passed to limit() */ + 0.0 /* Value that limit() is not to exceed */ + ); + + /* Duplicate in the CAM clip rspl if it exists */ + if (p->cclutTable != NULL) { + p->cclutTable->rev_set_limit( + p->cclutTable, /* this */ + p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0 ? icxLimitD_void : NULL, + /* Optional input space limit() function. */ + /* Function should evaluate in[0..di-1], and return */ + /* number that is not to exceed 0.0. NULL if not used. */ + (void *)p, /* Context passed to limit() */ + 0.0 /* Value that limit() is not to exceed */ + ); + } + + /* Figure Lmin and Lmax for icxKluma5 curve basis */ + if (setLminmax + && p->clutTable->di > p->clutTable->fdi) { /* If K generation makes sense */ + double wh[3], bk[3], kk[3]; + int mergeclut; /* Save/restore mergeclut value */ + + /* Get white/black in effective xlu PCS space */ + p->efv_wh_bk_points((icxLuBase *)p, wh, bk, kk); + + /* Convert from effective PCS (ie. Jab) to native XYZ or Lab PCS */ + mergeclut = p->mergeclut; /* Hack to be able to use inv_out_abs() */ + p->mergeclut = 0; /* if mergeclut is active. */ + icxLuLut_inv_out_abs(p, wh, wh); + icxLuLut_inv_out_abs(p, bk, bk); + icxLuLut_inv_out_abs(p, kk, kk); + p->mergeclut = mergeclut; /* Restore */ + + /* Convert to Lab PCS */ + if (p->natos == icSigXYZData) { /* Always do K rule in L space */ + icmXYZ2Lab(&icmD50, wh, wh); + icmXYZ2Lab(&icmD50, bk, bk); + icmXYZ2Lab(&icmD50, kk, kk); + } + p->Lmax = 0.01 * wh[0]; + if (p->ink.KonlyLmin != 0) + p->Lmin = 0.01 * kk[0]; + else + p->Lmin = 0.01 * bk[0]; + } else { + + /* Some sane defaults */ + p->Lmax = 1.0; + p->Lmin = 0.0; + } + + return 0; +} + + +/* Initialise the clut clipping information, ink limiting */ +/* and auxiliary parameter settings for all the rspl. */ +/* return 0 or error code */ +static int +setup_clip_icxLuLut( +icxLuLut *p /* Object being initialised */ +) { + double tmin[MXDIDO], tmax[MXDIDO]; + int i; + + /* Setup for inversion of multi-dim clut */ + + p->kch = -1; /* Not known yet */ + + /* Set auxiliaries */ + for (i = 0; i < p->inputChan; i++) + p->auxm[i] = 0; + + if (p->inputChan > p->outputChan) { + switch(p->natis) { + case icSigCmykData: + /* Should fix icm/xicc to remember K channel */ + p->auxm[3] = 1; /* K is the auxiliary channel */ + break; + default: + if (p->kch >= 0) /* It's been discovered */ + p->auxm[p->kch] = 1; + else { + p->pp->errc = 2; + sprintf(p->pp->err,"Unknown colorspace %s when setting auxliaries", + icm2str(icmColorSpaceSignature, p->natis)); + return p->pp->errc; + } + break; + } + } + +// p->auxlinf = NULL; /* Input space auxiliary linearisation function - not implemented */ +// p->auxlinf_ctx = NULL; /* Opaque context for auxliliary linearisation function */ + + /* Aproximate center of clut input gamut - used for */ + /* resolving multiple reverse solutions. */ + p->clutTable->get_in_range(p->clutTable, tmin, tmax); + for (i = 0; i < p->clutTable->di; i++) { + p->licent[i] = p->icent[i] = (tmin[i] + tmax[i])/2.0; + } + + /* Compute clip setup information relating to clut output gamut. */ + if (p->nearclip != 0 /* Near clip requested */ + || p->inputChan == 1) { /* or vector clip won't work */ + p->clip.nearclip = 1; + + } else { /* Vector clip */ + /* !!!! NOTE NOTE NOTE !!!! */ + /* We should re-write this to avoid calling p->clutTable->rev_interp(), */ + /* since this sets up all the rev acceleration tables for two calls, */ + /* if this lut is being setup from scattered data, and never used */ + /* for rev lookup */ + icColorSpaceSignature clutos = p->natos; + + fprintf(stderr, "!!!!! setup_clip_icxLuLut with vector clip - possibly unnecessary rev setup !!!!\n"); + p->clip.nearclip = 0; + p->clip.LabLike = 0; + p->clip.fdi = p->clutTable->fdi; + + switch(clutos) { + case icxSigJabData: + case icSigLabData: { + co pp; /* Room for all the solutions found */ + int nsoln; /* Number of solutions found */ + double cdir[MXDO]; /* Clip vector direction and length */ + + p->clip.LabLike = 1; + + /* Find high clip target */ + for (i = 0; i < p->inputChan; i++) + pp.p[i] = 0.0; /* Set aux values */ + pp.v[0] = 105.0; pp.v[1] = pp.v[2] = 0.0; /* PCS Target value */ + cdir[0] = cdir[1] = cdir[2] = 0.0; /* Clip Target */ + + p->inv_output(p, pp.v, pp.v); /* Compensate for output curve */ + p->inv_output(p, cdir, cdir); + + cdir[0] -= pp.v[0]; /* Clip vector */ + cdir[1] -= pp.v[1]; + cdir[2] -= pp.v[2]; + + /* PCS -> Device with clipping */ + nsoln = p->clutTable->rev_interp( + p->clutTable, /* rspl object */ + 0, /* No hint flags - might be in gamut, might vector clip */ + 1, /* Maxumum solutions to return */ + p->auxm, /* Auxiliary input targets */ + cdir, /* Clip vector direction and length */ + &pp); /* Input target and output solutions */ + /* returned solutions in pp[0..retval-1].p[] */ + nsoln &= RSPL_NOSOLNS; /* Get number of solutions */ + + if (nsoln != 1) { + p->pp->errc = 2; + sprintf(p->pp->err,"Failed to find high clip target for Lab space"); + return p->pp->errc; + } + + p->clip.ocent[0] = pp.v[0] - 0.001; /* Got main target */ + p->clip.ocent[1] = pp.v[1]; + p->clip.ocent[2] = pp.v[2]; + + /* Find low clip target */ + pp.v[0] = -5.0; pp.v[1] = pp.v[2] = 0.0; /* PCS Target value */ + cdir[0] = 100.0; cdir[1] = cdir[2] = 0.0; /* Clip Target */ + + p->inv_output(p, pp.v, pp.v); /* Compensate for output curve */ + p->inv_output(p, cdir, cdir); + cdir[0] -= pp.v[0]; /* Clip vector */ + cdir[1] -= pp.v[1]; + cdir[2] -= pp.v[2]; + + /* PCS -> Device with clipping */ + nsoln = p->clutTable->rev_interp( + p->clutTable, /* rspl object */ + RSPL_WILLCLIP, /* Since there was no locus, we expect to have to clip */ + 1, /* Maxumum solutions to return */ + NULL, /* No auxiliary input targets */ + cdir, /* Clip vector direction and length */ + &pp); /* Input target and output solutions */ + /* returned solutions in pp[0..retval-1].p[] */ + nsoln &= RSPL_NOSOLNS; /* Get number of solutions */ + if (nsoln != 1) { + p->pp->errc = 2; + sprintf(p->pp->err,"Failed to find low clip target for Lab space"); + return p->pp->errc; + } + + p->clip.ocentv[0] = pp.v[0] + 0.001 - p->clip.ocent[0]; /* Raw line vector */ + p->clip.ocentv[1] = pp.v[1] - p->clip.ocent[1]; + p->clip.ocentv[2] = pp.v[2] - p->clip.ocent[2]; + + /* Compute vectors length */ + for (p->clip.ocentl = 0.0, i = 0; i < 3; i++) + p->clip.ocentl += p->clip.ocentv[i] * p->clip.ocentv[i]; + p->clip.ocentl = sqrt(p->clip.ocentl); + if (p->clip.ocentl <= 1e-8) + p->clip.ocentl = 0.0; + + break; + } + case icSigXYZData: + // ~~~~~~1 need to add this. + + default: + /* Do a crude approximation, that may not work. */ + p->clutTable->get_out_range(p->clutTable, tmin, tmax); + for (i = 0; i < p->clutTable->fdi; i++) { + p->clip.ocent[i] = (tmin[i] + tmax[i])/2.0; + } + p->clip.ocentl = 0.0; + break; + } + } + return 0; +} + +/* Function to pass to rspl to set secondary input/output transfer functions */ +static void +icxLuLut_inout_func( + void *pp, /* icxLuLut */ + double *out, /* output value */ + double *in /* inut value */ +) { + icxLuLut *p = (icxLuLut *)pp; /* this */ + icmLuLut *luluto = (icmLuLut *)p->plu; /* Get icmLuLut object */ + double tin[MAX_CHAN]; + double tout[MAX_CHAN]; + int i; + + if (p->iol_out == 0) { /* fwd input */ +#ifdef INK_LIMIT_TEST + tout[p->iol_ch] = in[0]; +#else + for (i = 0; i < p->inputChan; i++) + tin[i] = 0.0; + tin[p->iol_ch] = in[0]; + luluto->input(luluto, tout, tin); +#endif + } else if (p->iol_out == 1) { /* fwd output */ + for (i = 0; i < p->outputChan; i++) + tin[i] = 0.0; + tin[p->iol_ch] = in[0]; + luluto->output(luluto, tout, tin); + } else { /* bwd input */ +#ifdef INK_LIMIT_TEST + tout[p->iol_ch] = in[0]; +#else + for (i = 0; i < p->inputChan; i++) + tin[i] = 0.0; + tin[p->iol_ch] = in[0]; + luluto->inv_input(luluto, tout, tin); + /* This won't be valid if matrix is used or there is a PCS conversion !!! */ + /* Shouldn't be a problem because this is only used for inverting dev->PCS ? */ + luluto->inv_in_abs(luluto, tout, tout); +#endif + } + out[0] = tout[p->iol_ch]; +} + +/* Function to pass to rspl to set clut up, when mergeclut is set */ +static void +icxLuLut_clut_merge_func( + void *pp, /* icxLuLut */ + double *out, /* output value */ + double *in /* input value */ +) { + icxLuLut *p = (icxLuLut *)pp; /* this */ + icmLuLut *luluto = (icmLuLut *)p->plu; /* Get icmLuLut object */ + + luluto->clut(luluto, out, in); + luluto->output(luluto, out, out); + luluto->out_abs(luluto, out, out); + + if (p->outs == icxSigJabData) { + p->cam->XYZ_to_cam(p->cam, out, out); + } +} + +/* Implimenation of Lut create from icc. */ +/* Note that xicc_get_luobj() will have set the pcsor & */ +/* intent to consistent values if Jab and/or icxAppearance */ +/* has been requested. */ +/* It will also have created the underlying icm lookup object */ +/* that is used to create and implement the icx one. The icm */ +/* will be used to translate from native to effective PCS, unless */ +/* the effective PCS is Jab, in which case the icm will be set to */ +/* have an effective PCS of XYZ. Since native<->effecive PCS conversion */ +/* is done at the to/from_abs() stage, none of it affects the individual */ +/* conversion steps, which will all talk the native PCS (unless merged). */ +static icxLuBase * +new_icxLuLut( +xicc *xicp, +int flags, /* clip, merge flags */ +icmLuBase *plu, /* Pointer to Lu we are expanding (ours) */ +icmLookupFunc func, /* Functionality requested */ +icRenderingIntent intent, /* Rendering intent */ +icColorSpaceSignature pcsor, /* PCS override (0 = def) */ +icxViewCond *vc, /* Viewing Condition (NULL if not using CAM) */ +icxInk *ink /* inking details (NULL for default) */ +) { + icxLuLut *p; /* Object being created */ + icmLuLut *luluto = (icmLuLut *)plu; /* Lookup Lut type object */ + + int i; + + /* Do basic creation and initialisation */ + if ((p = alloc_icxLuLut(xicp, plu, flags)) == NULL) + return NULL; + + p->func = func; + + /* Set LuLut "use" specific creation flags: */ + if (flags & ICX_CLIP_NEAREST) + p->nearclip = 1; + + if (flags & ICX_MERGE_CLUT) + p->mergeclut = 1; + + if (flags & ICX_FAST_SETUP) + p->fastsetup = 1; + + /* We're only implementing this under specific conditions. */ + if (flags & ICX_CAM_CLIP + && func == icmFwd + && !(p->mergeclut != 0 && pcsor == icxSigJabData)) /* Don't need camclip if merged Jab */ + p->camclip = 1; + + if (flags & ICX_INT_SEPARATE) { +fprintf(stderr,"~1 Internal optimised 4D separations not yet implemented!\n"); + p->intsep = 1; + } + + /* Init the CAM model if it will be used */ + if (pcsor == icxSigJabData || p->camclip) { + 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(p->cam, p->vc.Ev, p->vc.Wxyz, p->vc.La, p->vc.Yb, p->vc.Lv, p->vc.Yf, + p->vc.Fxyz, XICC_USE_HK); + } else { + p->cam = NULL; + } + + /* Remember the effective intent */ + p->intent = intent; + + /* Get the effective spaces of underlying icm */ + plu->spaces(plu, &p->ins, NULL, &p->outs, NULL, NULL, NULL, NULL, &p->pcs, NULL); + + /* Override with pcsor */ + /* We assume that any profile that has a CIE color as a "device" color */ + /* intends it to stay that way, and not be overridden. */ + if (pcsor == icxSigJabData) { + p->pcs = pcsor; + + if (xicp->pp->header->deviceClass == icSigAbstractClass) { + p->ins = pcsor; + p->outs = pcsor; + + } else if (xicp->pp->header->deviceClass != icSigLinkClass) { + if (func == icmBwd || func == icmGamut || func == icmPreview) + p->ins = pcsor; + if (func == icmFwd || func == icmPreview) + p->outs = pcsor; + } + } + + /* In general the native and effective ranges of the icx will be the same as the */ + /* underlying icm lookup object. */ + p->plu->get_lutranges(p->plu, p->ninmin, p->ninmax, p->noutmin, p->noutmax); + p->plu->get_ranges(p->plu, p->inmin, p->inmax, p->outmin, p->outmax); + + /* 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; + } + + /* If we have a merged clut, reflect this in the icx native PCS range. */ + /* Merging merges output processing (irrespective of whether we are using */ + /* the forward or backward cluts) */ + if (p->mergeclut != 0) { + int i; + for (i = 0; i < p->outputChan; i++) { + p->noutmin[i] = p->outmin[i]; + p->noutmax[i] = p->outmax[i]; + } + } + + /* ------------------------------- */ + /* Create rspl based input lookups */ + for (i = 0; i < p->inputChan; i++) { + if ((p->inputTable[i] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of input table rspl failed"); + p->del((icxLuBase *)p); + return NULL; + } + p->iol_out = 0; /* Input lookup */ + p->iol_ch = i; /* Chanel */ + p->inputTable[i]->set_rspl(p->inputTable[i], RSPL_NOFLAGS, + (void *)p, icxLuLut_inout_func, + &p->ninmin[i], &p->ninmax[i], (int *)&p->lut->inputEnt, &p->ninmin[i], &p->ninmax[i]); + } + + /* Setup center clip target for inversion */ + for (i = 0; i < p->inputChan; i++) { + p->inputClipc[i] = (p->ninmin[i] + p->ninmax[i])/2.0; + } + + /* Create rspl based reverse input lookups used in ink limit and locus range functions. */ + for (i = 0; i < p->inputChan; i++) { + int gres = p->inputTable[i]->g.mres; /* Same res as curve being inverted */ + if (gres < 256) + gres = 256; + if ((p->revinputTable[i] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of reverse input table rspl failed"); + p->del((icxLuBase *)p); + return NULL; + } + p->iol_out = 2; /* Input lookup */ + p->iol_ch = i; /* Chanel */ + p->revinputTable[i]->set_rspl(p->revinputTable[i], RSPL_NOFLAGS, + (void *)p, icxLuLut_inout_func, + &p->ninmin[i], &p->ninmax[i], &gres, &p->ninmin[i], &p->ninmax[i]); + } + + /* ------------------------------- */ + { + int gres[MXDI]; + + for (i = 0; i < p->inputChan; i++) + gres[i] = p->lut->clutPoints; + + /* Create rspl based multi-dim table */ + if ((p->clutTable = new_rspl((p->fastsetup ? RSPL_FASTREVSETUP : RSPL_NOFLAGS) + | (flags & ICX_VERBOSE ? RSPL_VERBOSE : RSPL_NOFLAGS), + p->inputChan, p->outputChan)) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of clut table rspl failed"); + p->del((icxLuBase *)p); + return NULL; + } + + if (p->mergeclut == 0) { /* Do this if it's not merged with clut, */ + p->clutTable->set_rspl(p->clutTable, RSPL_NOFLAGS, + (void *)luluto, (void (*)(void *, double *, double *))luluto->clut, + p->ninmin, p->ninmax, gres, p->noutmin, p->noutmax); + + } else { /* If mergeclut */ + p->clutTable->set_rspl(p->clutTable, RSPL_NOFLAGS, + (void *)p, icxLuLut_clut_merge_func, + p->ninmin, p->ninmax, gres, p->noutmin, p->noutmax); + + } + + /* clut clipping is setup separately */ + } + + /* ------------------------------- */ + /* Create rspl based output lookups */ + for (i = 0; i < p->outputChan; i++) { + if ((p->outputTable[i] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of output table rspl failed"); + p->del((icxLuBase *)p); + return NULL; + } + p->iol_out = 1; /* Output lookup */ + p->iol_ch = i; /* Chanel */ + p->outputTable[i]->set_rspl(p->outputTable[i], RSPL_NOFLAGS, + (void *)p, icxLuLut_inout_func, + &p->noutmin[i], &p->noutmax[i], (int *)&p->lut->outputEnt, &p->noutmin[i], &p->noutmax[i]); + } + + /* Setup center clip target for inversion */ + for (i = 0; i < p->outputChan; i++) { + p->outputClipc[i] = (p->noutmin[i] + p->noutmax[i])/2.0; + } + + /* ------------------------------- */ + + /* Setup all the clipping, ink limiting and auxiliary stuff, */ + /* in case a reverse call is used. Only do this if we know */ + /* the reverse stuff isn't going to fail due to channel limits. */ + if (p->clutTable->within_restrictedsize(p->clutTable)) { + + if (setup_ink_icxLuLut(p, ink, 1) != 0) { + p->del((icxLuBase *)p); + return NULL; + } + + if (setup_clip_icxLuLut(p) != 0) { + p->del((icxLuBase *)p); + return NULL; + } + } + + return (icxLuBase *)p; +} + + +/* Function to pass to rspl to set clut up, when camclip is going to be used. */ +/* We use the temporary icm fwd absolute xyz lookup, then convert to CAM Jab. */ +static void +icxLuLut_clut_camclip_func( + void *pp, /* icxLuLut */ + double *out, /* output value */ + double *in /* inut value */ +) { + icxLuLut *p = (icxLuLut *)pp; /* this */ + icmLuLut *luluto = (icmLuLut *)p->absxyzlu; + + luluto->clut(luluto, out, in); + luluto->output(luluto, out, out); + luluto->out_abs(luluto, out, out); + p->cam->XYZ_to_cam(p->cam, out, out); +} + +/* Initialise the additional CAM space clut rspl, used to compute */ +/* reverse lookup CAM clipping results when the camclip flag is set. */ +/* return error code. */ +static int +icxLuLut_init_clut_camclip( +icxLuLut *p) { + int e, gres[MXDI]; + + /* Setup so clut contains transform to CAM Jab */ + /* (camclip is only used in fwd or invfwd direction lookup) */ + double cmin[3], cmax[3]; + cmin[0] = 0.0; cmax[0] = 100.0; /* Nominal Jab output ranges */ + cmin[1] = -128.0; cmax[1] = 128.0; + cmin[2] = -128.0; cmax[2] = 128.0; + + /* Get icm lookup we need for seting up and using CAM icx clut */ + if ((p->absxyzlu = p->pp->pp->get_luobj(p->pp->pp, icmFwd, icAbsoluteColorimetric, + icSigXYZData, icmLuOrdNorm)) == NULL) { + p->pp->errc = p->pp->pp->errc; /* Copy error to xicc */ + strcpy(p->pp->err, p->pp->pp->err); + return p->pp->errc; + } + + /* Create CAM rspl based multi-dim table */ + if ((p->cclutTable = new_rspl((p->fastsetup ? RSPL_FASTREVSETUP : RSPL_NOFLAGS) + | (p->flags & ICX_VERBOSE ? RSPL_VERBOSE : RSPL_NOFLAGS), + p->inputChan, p->outputChan)) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of clut table rspl failed"); + return p->pp->errc; + } + + for (e = 0; e < p->inputChan; e++) + gres[e] = p->lut->clutPoints; + + /* Setup our special CAM space rspl */ + p->cclutTable->set_rspl(p->cclutTable, RSPL_NOFLAGS, (void *)p, + icxLuLut_clut_camclip_func, + p->ninmin, p->ninmax, gres, cmin, cmax); + + /* Duplicate the ink limit information for any reverse interpolation. */ + p->cclutTable->rev_set_limit( + p->cclutTable, /* this */ + p->ink.tlimit >= 0.0 || p->ink.klimit >= 0.0 ? icxLimitD_void : NULL, + /* Optional input space limit() function. */ + /* Function should evaluate in[0..di-1], and return */ + /* number that is not to exceed 0.0. NULL if not used. */ + (void *)p, /* Context passed to limit() */ + 0.0 /* Value that limit() is not to exceed */ + ); + return 0; +} + +/* ========================================================== */ +/* xicc creation code */ +/* ========================================================== */ + +/* Callback for computing delta E squared for two output (PCS) */ +/* values, passed as a callback to xfit */ +static double xfit_to_de2(void *cntx, double *in1, double *in2) { + icxLuLut *p = (icxLuLut *)cntx; + double rv; + + if (p->pcs == icSigLabData) { +#ifdef USE_CIE94_DE + rv = icmCIE94sq(in1, in2); +#else + rv = icmLabDEsq(in1, in2); +#endif + } else { + double lab1[3], lab2[3]; + icmXYZ2Lab(&icmD50, lab1, in1); + icmXYZ2Lab(&icmD50, lab2, in2); +#ifdef USE_CIE94_DE + rv = icmCIE94sq(lab1, lab2); +#else + rv = icmLabDEsq(lab1, lab2); +#endif + } + return rv; +} + +/* Same as above plus partial derivatives */ +static double xfit_to_dde2(void *cntx, double dout[2][MXDIDO], double *in1, double *in2) { + icxLuLut *p = (icxLuLut *)cntx; + double rv; + + if (p->pcs == icSigLabData) { + int j,k; + double tdout[2][3]; +#ifdef USE_CIE94_DE + rv = icxdCIE94sq(tdout, in1, in2); +#else + rv = icxdLabDEsq(tdout, in1, in2); +#endif + for (k = 0; k < 2; k++) { + for (j = 0; j < 3; j++) + dout[k][j] = tdout[k][j]; + } + } else { + double lab1[3], lab2[3]; + double dout12[2][3][3]; + double tdout[2][3]; + int i,j,k; + + icxdXYZ2Lab(&icmD50, lab1, dout12[0], in1); + icxdXYZ2Lab(&icmD50, lab2, dout12[1], in2); +#ifdef USE_CIE94_DE + rv = icxdCIE94sq(tdout, lab1, lab2); +#else + rv = icxdLabDEsq(tdout, lab1, lab2); +#endif + /* Compute partial derivative (is this correct ??) */ + for (k = 0; k < 2; k++) { + for (j = 0; j < 3; j++) { + dout[k][j] = 0.0; + for (i = 0; i < 3; i++) { + dout[k][j] += tdout[k][i] * dout12[k][i][j]; + } + } + } + } + return rv; +} + +#ifdef NEVER +/* Check partial derivative function within xfit_to_dde2() */ + +static double _xfit_to_dde2(void *cntx, double dout[2][MXDIDO], double *in1, double *in2) { + icxLuLut *pp = (icxLuLut *)cntx; + int k, i; + double rv, drv; + double trv; + + rv = xfit_to_de2(cntx, in1, in2); + drv = xfit_to_dde2(cntx, dout, in1, in2); + + if (fabs(rv - drv) > 1e-6) + printf("######## DDE2: RV MISMATCH is %f should be %f ########\n",rv,drv); + + /* Check each parameter delta */ + for (k = 0; k < 2; k++) { + for (i = 0; i < 3; i++) { + double *in; + double del; + + if (k == 0) + in = in1; + else + in = in2; + + in[i] += 1e-9; + trv = xfit_to_de2(cntx, in1, in2); + in[i] -= 1e-9; + + /* Check that del is correct */ + del = (trv - rv)/1e-9; + if (fabs(dout[k][i] - del) > 0.04) { + printf("######## DDE2: EXCESSIVE at in[%d][%d] is %f should be %f ########\n",k,i,dout[k][i],del); + } + } + } + return rv; +} + +#define xfit_to_dde2 _xfit_to_dde2 + +#endif + +/* Context for rspl setting input and output curves */ +typedef struct { + int iix; + int oix; + xfit *xf; /* Optimisation structure */ +} curvectx; + +/* Function to pass to rspl to set input and output */ +/* transfer function for xicc lookup function */ +static void +set_linfunc( + void *cc, /* curvectx structure */ + double *out, /* Device output value */ + double *in /* Device input value */ +) { + curvectx *c = (curvectx *)cc; /* this */ + xfit *p = c->xf; + + if (c->iix >= 0) { /* Input curve */ + *out = p->incurve(p, *in, c->iix); + } else if (c->oix >= 0) { /* Output curve */ + *out = p->outcurve(p, *in, c->oix); + } +} + +/* Function to pass to rspl to set inverse input transfer function, */ +/* used for ink limiting calculation. */ +static void +icxLuLut_invinput_func( + void *cc, /* curvectx structure */ + double *out, /* Device output value */ + double *in /* Device input value */ +) { + curvectx *c = (curvectx *)cc; /* this */ + xfit *p = c->xf; + + *out = p->invincurve(p, *in, c->iix); +} + + +/* Functions to pass to icc settables() to setup icc A2B Lut: */ + +/* Input table */ +static void set_input(void *cntx, double *out, double *in) { + icxLuLut *p = (icxLuLut *)cntx; + + if (p->noisluts != 0 && p->noipluts != 0) { /* Input table must be linear */ + int i; + for (i = 0; i < p->inputChan; i++) + out[i] = in[i]; + } else { + if (p->input(p, out, in) > 1) + error ("%d, %s",p->pp->errc,p->pp->err); + } +} + +/* clut */ +static void set_clut(void *cntx, double *out, double *in) { + icxLuLut *p = (icxLuLut *)cntx; + int f; + + if (p->clut(p, out, in) > 1) + error ("%d, %s",p->pp->errc,p->pp->err); + + /* Convert from efective pcs to natural pcs */ + if (p->pcs != p->plu->icp->header->pcs) { + if (p->pcs == icSigLabData) + icmLab2XYZ(&icmD50, out, out); + else + icmXYZ2Lab(&icmD50, out, out); + } +} + +/* output */ +static void set_output(void *cntx, double *out, double *in) { + icxLuLut *p = (icxLuLut *)cntx; + + if (p->nooluts != 0) { /* Output table must be linear */ + int i; + for (i = 0; i < p->outputChan; i++) + out[i] = in[i]; + } else { + if (p->output(p, out, in) > 1) + error ("%d, %s",p->pp->errc,p->pp->err); + } +} + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Routine to figure out a suitable black point for CMYK */ + +/* Structure to hold optimisation information */ +typedef struct { + icxLuLut *p; /* Object being created */ + double toAbs[3][3]; /* To abs from aprox relative */ + double p1[3]; /* white pivot point in abs Lab */ + double p2[3]; /* Point on vector towards black */ + double toll; /* Tollerance of black direction */ +} bfinds; + +/* Optimise device values to minimise L, while remaining */ +/* within the ink limit, and staying in line between p1 (white) and p2 (black dir) */ +static double bfindfunc(void *adata, double pv[]) { + bfinds *b = (bfinds *)adata; + double rv = 0.0; + double tt[3], Lab[3]; + co bcc; + double lr, ta, tb, terr; /* L ratio, target a, target b, target error */ + double ovr; + +//printf("~1 bfindfunc got %f %f %f %f\n", pv[0], pv[1], pv[2], pv[3]); + /* See if over ink limit or outside device range */ + ovr = icxLimit(b->p, pv); /* > 0.0 if outside device gamut or ink limit */ + if (ovr < 0.0) + ovr = 0.0; +//printf("~1 bfindfunc got ovr %f\n", ovr); + + /* Compute the absolute Lab value: */ + b->p->input(b->p, bcc.p, pv); /* Through input tables */ + b->p->clutTable->interp(b->p->clutTable, &bcc); /* Through clut */ + b->p->output(b->p, bcc.v, bcc.v); /* Through the output tables */ + + if (b->p->pcs != icSigXYZData) /* Convert PCS to XYZ */ + icmLab2XYZ(&icmD50, bcc.v, bcc.v); + + /* Convert from relative to Absolute colorimetric */ + icmMulBy3x3(tt, b->toAbs, bcc.v); + icmXYZ2Lab(&icmD50, Lab, tt); /* Convert to Lab */ + +#ifdef DEBUG +printf("~1 p1 = %f %f %f, p2 = %f %f %f\n",b->p1[0],b->p1[1],b->p1[2],b->p2[0],b->p2[1],b->p2[2]); +printf("~1 device value %f %f %f %f, Lab = %f %f %f\n",pv[0],pv[1],pv[2],pv[3],Lab[0],Lab[1],Lab[2]); +#endif + + /* Primary aim is to minimise L value */ + rv = Lab[0]; + + /* See how out of line from p1 to p2 we are */ + lr = (Lab[0] - b->p1[0])/(b->p2[0] - b->p1[0]); /* Distance towards p2 from p1 */ + ta = lr * (b->p2[1] - b->p1[1]) + b->p1[1]; /* Target a value */ + tb = lr * (b->p2[2] - b->p1[2]) + b->p1[2]; /* Target b value */ + + terr = (ta - Lab[1]) * (ta - Lab[1]) + + (tb - Lab[2]) * (tb - Lab[2]); + + if (terr < b->toll) /* Tollerance error doesn't count until it's over tollerance */ + terr = 0.0; + +#ifdef DEBUG +printf("~1 target error %f\n",terr); +#endif + rv += XICC_BLACK_FIND_ABERR_WEIGHT * terr; /* Make ab match more important than min. L */ + +#ifdef DEBUG +printf("~1 out of range error %f\n",ovr); +#endif + rv += 200 * ovr; + +#ifdef DEBUG +printf("~1 black find tc ret %f\n",rv); +#endif + return rv; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/* Create icxLuLut and underlying fwd Lut from scattered data */ +/* The scattered data is assumed to map Device -> native PCS */ +/* NOTE:- in theory once this icxLuLut is setup, it can be */ +/* called to translate color values. In practice I suspect */ +/* that the icxLuLut hasn't been setup completely enough to allows this. */ +/* Might be easier to close it and re-open it ? */ +static icxLuBase * +set_icxLuLut( +xicc *xicp, +icmLuBase *plu, /* Pointer to Lu we are expanding (ours) */ +icmLookupFunc func, /* Functionality requested */ +icRenderingIntent intent, /* Rendering intent */ +int flags, /* white/black point flags etc. */ +int nodp, /* Number of points */ +int nodpbw, /* Number of points to look for white & black patches in */ +cow *ipoints, /* Array of input points (Lab or XYZ normalized to 1.0) */ +icxMatrixModel *skm, /* Optional skeleton model (used for input profiles) */ +double dispLuminance, /* > 0.0 if display luminance value and is known */ +double wpscale, /* > 0.0 if white point is to be scaled */ +double smooth, /* RSPL smoothing factor, -ve if raw */ +double avgdev, /* reading Average Deviation as a prop. of the input range */ +icxViewCond *vc, /* Viewing Condition (NULL if not using CAM) */ +icxInk *ink, /* inking details (NULL for default) */ +int quality /* Quality metric, 0..3 */ +) { + icxLuLut *p; /* Object being created */ + icc *icco = xicp->pp; /* Underlying icc object */ + int luflags = 0; /* icxLuLut alloc clip, merge flags */ + int pcsy; /* Effective PCS L or Y chanel index */ + double pcsymax; /* Effective PCS L or Y maximum value */ + icmHeader *h = icco->header; /* Pointer to icc header */ + int maxchan; /* max(inputChan, outputChan) */ + int rsplflags = RSPL_NOFLAGS; /* Flags for scattered data rspl */ + int e, f, i, j; + double dwhite[MXDI], dblack[MXDI]; /* Device white and black values */ + double wp[3]; /* Absolute White point in XYZ */ + double bp[3]; /* Absolute Black point in XYZ */ + double dgwhite[MXDI]; /* Device space gamut boundary white (ie. RGB 1,1,1) */ + double oavgdev[MXDO]; /* Average output value deviation */ + int gres[MXDI]; /* RSPL/CLUT resolution */ + xfit *xf = NULL; /* Curve fitting class instance */ + + if (flags & ICX_VERBOSE) + rsplflags |= RSPL_VERBOSE; + +// if (flags & ICX_2PASSSMTH) +// rsplflags |= RSPL_2PASSSMTH; /* Smooth data using Gaussian */ + +// if (flags & ICX_EXTRA_FIT) +// rsplflags |= RSPL_EXTRAFIT2; /* Try extra hard to fit data */ + + luflags = flags; /* Transfer straight though ? */ + + /* Do basic creation and initialisation */ + if ((p = alloc_icxLuLut(xicp, plu, luflags)) == NULL) + return NULL; + + p->func = func; + + /* Set LuLut "use" specific creation flags: */ + if (flags & ICX_CLIP_NEAREST) + p->nearclip = 1; + + /* Set LuLut "create" specific flags: */ + if (flags & ICX_NO_IN_SHP_LUTS) + p->noisluts = 1; + + if (flags & ICX_NO_IN_POS_LUTS) + p->noipluts = 1; + + if (flags & ICX_NO_OUT_LUTS) + p->nooluts = 1; + + /* Get the effective spaces of underlying icm, and set icx the same */ + plu->spaces(plu, &p->ins, NULL, &p->outs, NULL, NULL, &p->intent, NULL, &p->pcs, NULL); + + /* For set_icx the effective pcs has to be the same as the native pcs */ + + if (p->pcs == icSigXYZData) { + pcsy = 1; /* Y of XYZ */ + pcsymax = 1.0; + } else { + pcsy = 0; /* L or Lab */ + pcsymax = 100.0; + } + + maxchan = p->inputChan > p->outputChan ? p->inputChan : p->outputChan; + + /* Translate overall average deviation into output channel deviation */ + /* (This is for rspl scattered data fitting smoothness adjustment) */ + /* (This could do with more tuning) */ + + /* For some unknown reason XYZ seems masively under-smoothed, so bump it up */ + if (p->pcs == icSigXYZData) { + oavgdev[0] = 1.0 * 0.60 * avgdev; + oavgdev[1] = 1.0 * 1.00 * avgdev; + oavgdev[2] = 1.0 * 0.60 * avgdev; + } else if (p->pcs == icSigLabData) { + oavgdev[0] = 1.00 * avgdev; + oavgdev[1] = 0.60 * avgdev; + oavgdev[2] = 0.60 * avgdev; + } else { /* Hmmm */ + for (f = 0; f < p->outputChan; f++) + oavgdev[f] = avgdev; + } + + /* In general the native and effective ranges of the icx will be the same as the */ + /* underlying icm lookup object. */ + p->plu->get_lutranges(p->plu, p->ninmin, p->ninmax, p->noutmin, p->noutmax); + p->plu->get_ranges(p->plu, p->inmin, p->inmax, p->outmin, p->outmax); + + /* ??? Does this ever happen with set_icxLuLut() ??? */ + /* 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; + } + + /* ------------------------------- */ + + if (flags & ICX_VERBOSE) + printf("Estimating white point\n"); + + icmXYZ2Ary(wp, icmD50); /* Set a default value - D50 */ + icmXYZ2Ary(bp, icmBlack); /* Set a default value - absolute black */ + + { + /* Figure out as best we can the device white and black points */ + + /* Compute device white and black points as if */ + /* we are doing an Output or Display device */ + { + switch (h->colorSpace) { + + case icSigCmykData: + for (e = 0; e < p->inputChan; e++) { + dwhite[e] = 0.0; + dblack[e] = 1.0; + } + break; + case icSigCmyData: + for (e = 0; e < p->inputChan; e++) { + dwhite[e] = 0.0; + dblack[e] = 1.0; + } + break; + case icSigRgbData: + for (e = 0; e < p->inputChan; e++) { + dwhite[e] = 1.0; + dblack[e] = 0.0; + } + break; + + case icSigGrayData: { /* Could be additive or subtractive */ + double maxY, minY; /* Y min and max */ + double maxd, mind; /* Corresponding device values */ + + maxY = -1e8, minY = 1e8; + + /* Figure out if it's additive or subtractive */ + for (i = 0; i < nodpbw; i++) { + double Y; + if (p->pcs != icSigXYZData) { /* Convert white point to XYZ */ + double xyz[3]; + icmLab2XYZ(&icmD50, xyz, ipoints[i].v); + Y = xyz[1]; + } else { + Y = ipoints[i].v[1]; + } + + if (Y > maxY) { + maxY = Y; + maxd = ipoints[i].p[0]; + } + if (Y < minY) { + minY = Y; + mind = ipoints[i].p[0]; + } + } + if (maxd < mind) { /* Subtractive */ + for (e = 0; e < p->inputChan; e++) { + dwhite[e] = 0.0; + dblack[e] = 1.0; + } + } else { /* Additive */ + for (e = 0; e < p->inputChan; e++) { + dwhite[e] = 1.0; + dblack[e] = 0.0; + } + } + break; + } + default: + xicp->errc = 1; + sprintf(xicp->err,"set_icxLuLut: can't handle color space %s", + icm2str(icmColorSpaceSignature, h->colorSpace)); + p->del((icxLuBase *)p); + return NULL; + break; + } + } + /* dwhite is what we want for dgwhite[], used for XFIT_OUT_WP_REL_US */ + for (e = 0; e < p->inputChan; e++) + dgwhite[e] = dwhite[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 */ + if (p->pcs == icSigXYZData) + icmXYZ2Lab(&icmD50, labv, ipoints[i].v); + else + icmCpy3(labv,ipoints[i].v); + +#ifdef NEVER + /* Choose largest Y or L* */ + if (ipoints[i].v[pcsy] > wpy) { + wp[0] = ipoints[i].v[0]; + wp[1] = ipoints[i].v[1]; + wp[2] = ipoints[i].v[2]; + for (e = 0; e < p->inputChan; e++) + dwhite[e] = ipoints[i].p[e]; + wpy = ipoints[i].v[pcsy]; + wix = i; + } +#else + + /* Tilt things towards D50 neutral white patches */ + yv = labv[0] - 0.3 * sqrt(labv[1] * labv[1] + labv[2] * labv[2]); +//printf("~1 patch %d Lab = %s, yv = %f\n",i+1,icmPdv(3,labv),yv); + if (yv > wpy) { +//printf("~1 best so far\n"); + wp[0] = ipoints[i].v[0]; + wp[1] = ipoints[i].v[1]; + wp[2] = ipoints[i].v[2]; + for (e = 0; e < p->inputChan; e++) + dwhite[e] = ipoints[i].p[e]; + wpy = yv; + wix = i; + } +#endif + if (ipoints[i].v[pcsy] < bpy) { + bp[0] = ipoints[i].v[0]; + bp[1] = ipoints[i].v[1]; + bp[2] = ipoints[i].v[2]; + for (e = 0; e < p->inputChan; e++) + dblack[e] = ipoints[i].p[e]; + bpy = ipoints[i].v[pcsy]; + bix = i; + } + } + /* Make sure values are XYZ */ + if (p->pcs != icSigXYZData) { + icmLab2XYZ(&icmD50, wp, wp); + icmLab2XYZ(&icmD50, bp, bp); + } + + if (flags & ICX_VERBOSE) { + printf("Picked white patch %d with dev = %s\n XYZ = %s, Lab = %s\n", + wix+1, icmPdv(p->inputChan, dwhite), icmPdv(3, wp), icmPLab(wp)); + printf("Picked black patch %d with dev = %s\n XYZ = %s, Lab = %s\n", + bix+1, icmPdv(p->inputChan, dblack), icmPdv(3, bp), icmPLab(bp)); + } + + /* else Output or Display device */ + } else { + /* We assume that the output 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 icSigCmykData: + 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 + && ipoints[i].p[3] < 0.001) { + wp[0] += ipoints[i].v[0]; + wp[1] += ipoints[i].v[1]; + wp[2] += ipoints[i].v[2]; + nw++; + } + } + break; + 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; + + case icSigGrayData: { /* Could be additive or subtractive */ + double minwp[3], maxwp[3]; + int nminwp = 0, nmaxwp = 0; + + minwp[0] = minwp[1] = minwp[2] = 0.0; + maxwp[0] = maxwp[1] = maxwp[2] = 0.0; + + /* Look for both */ + for (i = 0; i < nodpbw; i++) { + if (ipoints[i].p[0] < 0.001) + minwp[0] += ipoints[i].v[0]; + minwp[1] += ipoints[i].v[1]; + minwp[2] += ipoints[i].v[2]; { + nminwp++; + } + if (ipoints[i].p[0] > 0.999) + maxwp[0] += ipoints[i].v[0]; + maxwp[1] += ipoints[i].v[1]; + maxwp[2] += ipoints[i].v[2]; { + nmaxwp++; + } + } + if (nminwp > 0) { /* Subtractive */ + wp[0] = minwp[0]; + wp[1] = minwp[1]; + wp[2] = minwp[2]; + nw = nminwp; + if (minwp[pcsy]/nminwp < (0.5 * pcsymax)) + nw = 0; /* Looks like a mistake */ + } + if (nmaxwp > 0 /* Additive */ + && (nminwp == 0 || maxwp[pcsy]/nmaxwp > minwp[pcsy]/nminwp)) { + wp[0] = maxwp[0]; + wp[1] = maxwp[1]; + wp[2] = maxwp[2]; + nw = nmaxwp; + if (maxwp[pcsy]/nmaxwp < (0.5 * pcsymax)) + nw = 0; /* Looks like a mistake */ + } + break; + } + + default: + xicp->errc = 1; + sprintf(xicp->err,"set_icxLuLut: can't handle color space %s", + icm2str(icmColorSpaceSignature, h->colorSpace)); + p->del((icxLuBase *)p); + return NULL; + break; + } + + if (nw == 0) { + xicp->errc = 1; + sprintf(xicp->err,"set_icxLuLut: can't handle test points without a white patch"); + p->del((icxLuBase *)p); + return NULL; + } + wp[0] /= (double)nw; + wp[1] /= (double)nw; + wp[2] /= (double)nw; + if (p->pcs != icSigXYZData) /* Convert white point to XYZ */ + icmLab2XYZ(&icmD50, wp, wp); + } + + if (flags & ICX_VERBOSE) { + printf("Approximate White point XYZ = %s, Lab = %s\n", icmPdv(3,wp),icmPLab(wp)); + } + } + + if (h->colorSpace == icSigGrayData) { /* Don't use device or PCS curves for monochrome */ + p->noisluts = p->noipluts = p->nooluts = 1; + } + + if ((flags & ICX_VERBOSE) && (p->noisluts == 0 || p->noipluts == 0 || p->nooluts == 0)) + printf("Creating optimised per channel curves\n"); + + /* Set the target CLUT grid resolution so in/out curves can be optimised for it */ + for (e = 0; e < p->inputChan; e++) + gres[e] = p->lut->clutPoints; + + /* Setup and then create xfit object that does most of the work */ + { + int xfflags = 0; /* xfit flags */ + double in_min[MXDI]; /* Input value scaling minimum */ + double in_max[MXDI]; /* Input value scaling maximum */ + double out_min[MXDO]; /* Output value scaling minimum */ + double out_max[MXDO]; /* Output value scaling maximum */ + int iluord, sluord, oluord; + int iord[MXDI]; /* Input curve orders */ + int sord[MXDI]; /* Input sub-grid curve orders */ + int oord[MXDO]; /* Output curve orders */ + double shp_smooth[MXDI];/* Smoothing factors for each curve, nom = 1.0 */ + double out_smooth[MXDO]; + + optcomb tcomb = oc_ipo; /* Create all by default */ + + if ((xf = CAT2(new_, xfit)()) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of xfit object failed"); + p->del((icxLuBase *)p); + return NULL; + } + + /* Setup for optimising run */ + if (p->noisluts) + tcomb &= ~oc_i; + + if (p->noipluts) + tcomb &= ~oc_p; + + if (p->nooluts) + tcomb &= ~oc_o; + + if (flags & ICX_VERBOSE) + xfflags |= XFIT_VERB; + + if (flags & ICX_SET_WHITE) { + + xfflags |= XFIT_OUT_WP_REL; + if ((flags & ICX_SET_WHITE_C) == ICX_SET_WHITE_C) { + xfflags |= XFIT_OUT_WP_REL_C; + } + else if ((flags & ICX_SET_WHITE_US) == ICX_SET_WHITE_US) + xfflags |= XFIT_OUT_WP_REL_US; + + if (p->pcs != icSigXYZData) + xfflags |= XFIT_OUT_LAB; + } + + /* If asked, clip the absolute white point to have Y <= 1.0 ? */ + if (flags & ICX_CLIP_WB) + xfflags |= XFIT_CLIP_WP; + + /* With current B2A code, make sure a & b curves */ + /* pass through zero. */ + if (p->pcs == icSigLabData) { + xfflags |= XFIT_OUT_ZERO; + } + + /* Let xfit create the clut */ + xfflags |= XFIT_MAKE_CLUT; + + /* Set the curve orders for input (device) and output (PCS) */ + if (quality >= 3) { /* Ultra high */ + iluord = 25; + sluord = 4; + oluord = 25; + } else if (quality == 2) { /* High */ + iluord = 20; + sluord = 3; + oluord = 20; + } else if (quality == 1) { /* Medium */ + iluord = 17; + sluord = 2; + oluord = 17; + } else { /* Low */ + iluord = 10; + sluord = 1; + oluord = 10; + } + for (e = 0; e < p->inputChan; e++) { + iord[e] = iluord; + sord[e] = sluord; + in_min[e] = p->inmin[e]; + in_max[e] = p->inmax[e]; + shp_smooth[e] = SHP_SMOOTH; + } + + for (f = 0; f < p->outputChan; f++) { + oord[f] = oluord; + out_min[f] = p->outmin[f]; + out_max[f] = p->outmax[f]; + + /* Hack to prevent a convex L curve pushing */ + /* the clut L values above the maximum value */ + /* that can be represented, causing clipping. */ + /* Do this by making sure that the L curve pivots */ + /* through 100.0 to 100.0 */ + if (f == 0 && p->pcs == icSigLabData) { + if (out_min[f] < 0.0001 && out_max[f] > 100.0) { + out_max[f] = 100.0; + } + } + out_smooth[f] = OUT_SMOOTH1; + if (f != 0 && p->pcs == icSigLabData) /* a* & b* */ + out_smooth[f] = OUT_SMOOTH2; + } + +//printf("~1 wp before xfit %f %f %f\n", wp[0], wp[1], wp[2]); + /* Create input, sub and output per channel curves (if configured), */ + /* adjust for white point to make relative (if configured), */ + /* and create clut rspl using xfit class. */ + /* The true white point for the returned curves and rspl is returned. */ + if (xf->fit(xf, xfflags, p->inputChan, p->outputChan, + rsplflags, wp, dwhite, wpscale, dgwhite, + ipoints, nodp, skm, in_min, in_max, gres, out_min, out_max, + smooth, oavgdev, iord, sord, oord, shp_smooth, out_smooth, tcomb, + (void *)p, xfit_to_de2, xfit_to_dde2) != 0) { + p->pp->errc = 2; + sprintf(p->pp->err,"xfit fitting failed"); + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } +//printf("~1 wp after xfit %f %f %f\n", wp[0], wp[1], wp[2]); + + /* - - - - - - - - - - - - - - - */ + /* Set the xicc input curve rspl */ + for (e = 0; e < p->inputChan; e++) { + curvectx cx; + + cx.xf = xf; + cx.oix = -1; + cx.iix = e; + + if ((p->inputTable[e] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of input table rspl failed"); + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + + p->inputTable[e]->set_rspl(p->inputTable[e], RSPL_NOFLAGS, + (void *)&cx, set_linfunc, + &p->ninmin[e], &p->ninmax[e], + (int *)&p->lut->inputEnt, + &p->ninmin[e], &p->ninmax[e]); + } + + /* - - - - - - - - - - - - - - - */ + /* Set the xicc output curve rspl */ + for (f = 0; f < p->outputChan; f++) { + curvectx cx; + + cx.xf = xf; + cx.iix = -1; + cx.oix = f; + + if ((p->outputTable[f] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of output table rspl failed"); + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + + p->outputTable[f]->set_rspl(p->outputTable[f], RSPL_NOFLAGS, + (void *)&cx, set_linfunc, + &p->noutmin[f], &p->noutmax[f], + (int *)&p->lut->outputEnt, + &p->noutmin[f], &p->noutmax[f]); + + } + } + + if (flags & ICX_VERBOSE) + printf("Creating fast inverse input lookups\n"); + + /* Create rspl based reverse input lookups used in ink limit function. */ + for (e = 0; e < p->inputChan; e++) { + int res = 256; + curvectx cx; + + cx.xf = xf; + cx.oix = -1; + cx.iix = e; + + if ((p->revinputTable[e] = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) { + p->pp->errc = 2; + sprintf(p->pp->err,"Creation of reverse input table rspl failed"); + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + p->revinputTable[e]->set_rspl(p->revinputTable[e], RSPL_NOFLAGS, + (void *)&cx, icxLuLut_invinput_func, + &p->ninmin[e], &p->ninmax[e], &res, &p->ninmin[e], &p->ninmax[e]); + } + + + /* ------------------------------- */ + /* Set clut lookup table from xfit */ + p->clutTable = xf->clut; + xf->clut = NULL; + + /* Setup all the clipping, ink limiting and auxiliary stuff, */ + /* in case a reverse call is used. Need to avoid relying on inking */ + /* stuff that makes use of the white/black points, since they haven't */ + /* been set up properly yet. */ + if (setup_ink_icxLuLut(p, ink, 0) != 0) { + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + + /* Deal with finalizing white/black points */ + { + + if ((flags & ICX_SET_WHITE) && (flags & ICX_VERBOSE)) { + printf("White point XYZ = %s, Lab = %s\n", icmPdv(3,wp),icmPLab(wp)); + } + + /* Lookup the black point */ + { /* Black Point Tag: */ + co bcc; + + if (flags & ICX_VERBOSE) + printf("Find black point\n"); + + /* !!! Hmm. For CMY and RGB we are simply using the device */ + /* combination values as the black point. In reality we might */ + /* want to have the option of using a neutral black point, */ + /* just like CMYK ?? */ + + /* For CMYK devices, we choose a black that has minumum L within */ + /* the ink limits, and if XICC_NEUTRAL_CMYK_BLACK it will in the direction */ + /* that has the same chromaticity as the white point, else choose the same */ + /* Lab vector direction as K, with the minimum L value. */ + /* (Note this is duplicated in xicc.c icxLu_comp_bk_point() !!!) */ + if (h->deviceClass != icSigInputClass + && h->colorSpace == icSigCmykData) { + bfinds bfs; /* Callback context */ + double sr[MXDO]; /* search radius */ + double tt[MXDO]; /* Temporary */ + double rs0[MXDO], rs1[MXDO]; /* Random start candidates */ + int trial; + double brv; + int kch = 3; + + /* Setup callback function context */ + bfs.p = p; + bfs.toll = XICC_BLACK_POINT_TOLL; + + /* !!! we should use an accessor funcion of xfit !!! */ + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + bfs.toAbs[i][j] = xf->toAbs[i][j]; + } + } + + /* Lookup abs Lab value of white point */ + icmXYZ2Lab(&icmD50, bfs.p1, wp); + +#ifdef XICC_NEUTRAL_CMYK_BLACK + icmScale3(tt, wp, 0.02); /* Scale white XYZ towards 0,0,0 */ + icmXYZ2Lab(&icmD50, bfs.p2, tt); /* Convert black direction to Lab */ + if (flags & ICX_VERBOSE) + printf("Neutral black direction (Lab) = %f %f %f\n",bfs.p2[0], bfs.p2[1], bfs.p2[2]); + +#else /* K direction black */ + /* Now figure abs Lab value of K only, as the direction */ + /* to use for the rich black. */ + for (e = 0; e < p->inputChan; e++) + bcc.p[e] = 0.0; + if (p->ink.klimit < 0.0) + bcc.p[kch] = 1.0; + else + bcc.p[kch] = p->ink.klimit; /* K value */ + + p->input(p, bcc.p, bcc.p); /* Through input tables */ + p->clutTable->interp(p->clutTable, &bcc); /* Through clut */ + p->output(p, bcc.v, bcc.v); /* Through the output tables */ + + if (p->pcs != icSigXYZData) /* Convert PCS to XYZ */ + icmLab2XYZ(&icmD50, bcc.v, bcc.v); + + /* Convert from relative to Absolute colorimetric */ + icmMulBy3x3(tt, xf->toAbs, bcc.v); + icmXYZ2Lab(&icmD50, bfs.p2, tt); /* Convert K only black point to Lab */ + + if (flags & ICX_VERBOSE) + printf("K only black direction (Lab) = %f %f %f\n",bfs.p2[0], bfs.p2[1], bfs.p2[2]); +#endif + /* Start with the K only as the current best value */ + brv = bfindfunc((void *)&bfs, bcc.p); +//printf("~1 initial brv for K only = %f\n",brv); + + /* Set the random start 0 location as 000K */ + /* and the random start 1 location as CMY0 */ + { + double tv; + + for (e = 0; e < p->inputChan; e++) + rs0[e] = 0.0; + if (p->ink.klimit < 0.0) + rs0[kch] = 1.0; + else + rs0[kch] = p->ink.klimit; /* K value */ + + if (p->ink.tlimit < 0.0) + tv = 1.0; + else + tv = p->ink.tlimit/(p->inputChan - 1.0); + for (e = 0; e < p->inputChan; e++) + rs1[e] = tv; + rs1[kch] = 0.0; /* K value */ + } + + /* Find the device black point using optimization */ + /* Do several trials to avoid local minima. */ + rand32(0x12345678); /* Make trial values deterministic */ + for (trial = 0; trial < 200; trial++) { + double rv; /* Temporary */ + + /* Start first trial at 000K */ + if (trial == 0) { + for (j = 0; j < p->inputChan; j++) { + tt[j] = rs0[j]; + sr[j] = 0.1; + } + + } else { + /* Base is random between 000K and CMY0: */ + if (trial < 100) { + rv = d_rand(0.0, 1.0); + for (j = 0; j < p->inputChan; j++) { + tt[j] = rv * rs0[j] + (1.0 - rv) * rs1[j]; + sr[j] = 0.1; + } + /* Base on current best */ + } else { + for (j = 0; j < p->inputChan; j++) { + tt[j] = bcc.p[j]; + sr[j] = 0.1; + } + } + + /* Then add random start offset */ + for (rv = 0.0, j = 0; j < p->inputChan; j++) { + tt[j] += d_rand(-0.5, 0.5); + if (tt[j] < 0.0) + tt[j] = 0.0; + else if (tt[j] > 1.0) + tt[j] = 1.0; + } + } + + /* Clip to be within ink limit */ + icxDoLimit(p, tt, tt); + + if (powell(&rv, p->inputChan, tt, sr, 0.000001, 1000, + bfindfunc, (void *)&bfs, NULL, NULL) == 0) { +//printf("~1 trial %d, rv %f bp %f %f %f %f\n",trial,rv,tt[0],tt[1],tt[2],tt[3]); + if (rv < brv) { +//printf("~1 new best\n"); + brv = rv; + for (j = 0; j < p->inputChan; j++) + bcc.p[j] = tt[j]; + } + } + } + if (brv > 1000.0) + error ("set_icxLuLut: Black point powell failed"); + + /* Make sure resulting device values are strictly in range */ + for (j = 0; j < p->inputChan; j++) { + if (bcc.p[j] < 0.0) + bcc.p[j] = 0.0; + else if (bcc.p[j] > 1.0) + bcc.p[j] = 1.0; + } + /* Now have device black in bcc.p[] */ +//printf("~1 Black point is CMYK %f %f %f %f\n", bcc.p[0], bcc.p[1], bcc.p[2], bcc.p[3]); + + /* Else not a CMYK output device, */ + /* use the previously determined device black value. */ + } else { + for (e = 0; e < p->inputChan; e++) + bcc.p[e] = dblack[e]; + } + + /* Lookup the PCS for the device black: */ + p->input(p, bcc.p, bcc.p); /* Through input tables */ + p->clutTable->interp(p->clutTable, &bcc); /* Through clut */ + p->output(p, bcc.v, bcc.v); /* Through the output tables */ + + if (p->pcs != icSigXYZData) /* Convert PCS to XYZ */ + icmLab2XYZ(&icmD50, bcc.v, bcc.v); + + /* Convert from relative to Absolute colorimetric */ + icmMulBy3x3(bp, xf->toAbs, bcc.v); + + /* Got XYZ black point in bp[] */ + if (flags & ICX_VERBOSE) { + printf("Black point XYZ = %s, Lab = %s\n", icmPdv(3,bp),icmPLab(bp)); + } + + /* Some ICC sw gets upset if the bp is at all -ve. */ + /* Clip it if this is the case */ + /* (Or we could get xfit to rescale the rspl instead ??) */ + if ((flags & ICX_CLIP_WB) + && (bp[0] < 0.0 || bp[1] < 0.0 || bp[2] < 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 this is a display, adjust the white point to be */ + /* exactly Y = 1.0, and compensate the dispLuminance */ + /* and black point accordingly. The Lut is already set to */ + /* assume device white maps to perfect PCS white. */ + if (h->deviceClass == icSigDisplayClass) { + double scale = 1.0/wp[1]; + int i; + + dispLuminance /= scale; + + for (i = 0; i < 3; i++) { + wp[i] *= scale; + bp[i] *= scale; + } + } + + if (h->deviceClass == icSigDisplayClass + && dispLuminance > 0.0) { + icmXYZArray *wo; + if ((wo = (icmXYZArray *)icco->read_tag( + icco, icSigLuminanceTag)) == NULL) { + xicp->errc = 1; + sprintf(xicp->err,"icx_set_luminance: couldn't find luminance tag"); + p->del((icxLuBase *)p); + return NULL; + } + if (wo->ttype != icSigXYZArrayType) { + xicp->errc = 1; + sprintf(xicp->err,"luminance: tag has wrong type"); + p->del((icxLuBase *)p); + return NULL; + } + + wo->size = 1; + wo->allocate((icmBase *)wo); /* Allocate space */ + wo->data[0].X = 0.0; + wo->data[0].Y = dispLuminance * wp[1]; + wo->data[0].Z = 0.0; + + if (flags & ICX_VERBOSE) + printf("Display Luminance = %f\n", wo->data[0].Y); + } + + /* Write white and black points */ + if (flags & ICX_SET_WHITE) { /* White Point Tag: */ + icmXYZArray *wo; + if ((wo = (icmXYZArray *)icco->read_tag( + icco, icSigMediaWhitePointTag)) == NULL) { + xicp->errc = 1; + sprintf(xicp->err,"icx_set_white_black: couldn't find white tag"); + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + if (wo->ttype != icSigXYZArrayType) { + xicp->errc = 1; + sprintf(xicp->err,"icx_set_white_black: white tag has wrong type"); + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + + wo->size = 1; + wo->allocate((icmBase *)wo); /* Allocate space */ + wo->data[0].X = wp[0]; + wo->data[0].Y = wp[1]; + wo->data[0].Z = wp[2]; + } + if (flags & ICX_SET_BLACK) { /* Black Point Tag: */ + icmXYZArray *wo; + if ((wo = (icmXYZArray *)icco->read_tag( + icco, icSigMediaBlackPointTag)) == NULL) { + xicp->errc = 1; + sprintf(xicp->err,"icx_set_white_black: couldn't find black tag"); + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + if (wo->ttype != icSigXYZArrayType) { + xicp->errc = 1; + sprintf(xicp->err,"icx_set_white_black: black tag has wrong type"); + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + + wo->size = 1; + wo->allocate((icmBase *)wo); /* Allocate space */ + wo->data[0].X = bp[0]; + wo->data[0].Y = bp[1]; + wo->data[0].Z = bp[2]; + } + if ((flags & ICX_SET_WHITE) || (flags & ICX_SET_BLACK)) { + /* Make sure ICC white/black point lookup notices the new white and black points */ + p->plu->init_wh_bk(p->plu); + } + + /* Setup the clut clipping, ink limiting and auxiliary stuff again */ + /* since re_set_rspl will have invalidated */ + if (setup_ink_icxLuLut(p, ink, 1) != 0) { + xf->del(xf); + p->del((icxLuBase *)p); + return NULL; + } + } + + /* Done with xfit now */ + xf->del(xf); + xf = NULL; + + if (setup_clip_icxLuLut(p) != 0) { + p->del((icxLuBase *)p); + return NULL; + } + + /* ------------------------------- */ + +//Debugging clipping of clut +//printf("~1 xlut.c: noutmin = %f %f %f\n", p->noutmin[0], p->noutmin[1], p->noutmin[2]); +//printf("~1 xlut.c: noutmax = %f %f %f\n", p->noutmax[0], p->noutmax[1], p->noutmax[2]); +//printf("~1 xlut.c: outmin = %f %f %f\n", p->outmin[0], p->outmin[1], p->outmin[2]); +//printf("~1 xlut.c: outmax = %f %f %f\n", p->outmax[0], p->outmax[1], p->outmax[2]); + + /* Do a specific test for out of Lab encoding range RGB primaries */ + /* (A more general check seems to get false positives - why ??) */ + if (h->pcs == icSigLabData + && ( h->deviceClass == icSigDisplayClass + || h->deviceClass == icSigOutputClass) + && h->colorSpace == icSigRgbData) { + double dev[3] = { 0.0, 0.0, 0.0 }; + double pcs[3]; + double clip = 0; + + for (i = 0; i < 3; i++) { + dev[i] = 1.0; + + if (p->clut(p, pcs, dev) > 1) + error ("%d, %s",p->pp->errc,p->pp->err); + + /* Convert from efective pcs to natural pcs */ + if (p->pcs != icSigLabData) + icmXYZ2Lab(&icmD50, pcs, pcs); + + if (pcs[1] < -128.0 || pcs[1] > 128.0 + || pcs[2] < -128.0 || pcs[2] > 128.0) { + warning("\n *** %s primary value can't be encoded in L*a*b* PCS (%f %f %f)", + i == 0 ? "Red" : i == 1 ? "Green" : "Blue", pcs[0],pcs[1],pcs[2]); + clip = 1; + } + dev[i] = 0.0; + } + if (clip) + a1logw(g_log," *** Try switching to XYZ PCS ***\n"); + } + + /* Use our rspl's to set the icc Lut AtoB table values. */ + /* Use helper function to do the hard work. */ + if (p->lut->set_tables(p->lut, ICM_CLUT_SET_EXACT, (void *)p, + h->colorSpace, /* Input color space */ + h->pcs, /* Output color space */ + set_input, /* Input transfer function, Dev->Dev' */ + NULL, NULL, /* Use default Maximum range of Dev' values */ + set_clut, /* Dev' -> PCS' transfer function */ + NULL, NULL, /* Use default Maximum range of PCS' values */ + set_output /* Linear output transform PCS'->PCS */ + ) != 0) + error("Setting 16 bit %s->%s Lut failed: %d, %s", + icm2str(icmColorSpaceSignature, h->colorSpace), + icm2str(icmColorSpaceSignature, h->pcs), + p->pp->pp->errc,p->pp->pp->err); + +#ifdef WARN_CLUT_CLIPPING + if (p->pp->pp->warnc) { + warning("Values clipped in setting A2B LUT!"); + if (p->pp->pp->warnc == 2 && h->pcs == icSigLabData) { + a1logw(g_log,"*** Try switching to XYZ PCS ***\n"); + } + } +#endif /* WARN_CLUT_CLIPPING */ + + /* Init a CAM model in case it will be used (ie. in profile with camclip flag) */ + 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(p->cam, p->vc.Ev, p->vc.Wxyz, p->vc.La, p->vc.Yb, p->vc.Lv, p->vc.Yf, + p->vc.Fxyz, XICC_USE_HK); + + if (flags & ICX_VERBOSE) + printf("Done A to B table creation\n"); + + return (icxLuBase *)p; +} + +/* ========================================================== */ +/* Gamut boundary code. */ +/* ========================================================== */ + + +/* Context for creating gamut boundary points fro, xicc */ +typedef struct { + gamut *g; /* Gamut being created */ + icxLuLut *x; /* xLut we are working from */ + icxLuBase *flu; /* Forward xlookup */ + double in[MAX_CHAN]; /* Device input values */ +} lutgamctx; + +/* Function to hand to zbrent to find a clut input' value at the ink limit */ +/* Returns value < 0.0 when within gamut, > 0.0 when out of gamut */ +static double icxLimitFind(void *fdata, double tp) { + int i; + double in[MAX_CHAN]; + lutgamctx *p = (lutgamctx *)fdata; + double tt; + + for (i = 0; i < p->x->inputChan; i++) { + in[i] = tp * p->in[i]; /* Scale given input value */ + } + + tt = icxLimitD(p->x, in); + + return tt; +} + +/* Function to pass to rspl to create gamut boundary from */ +/* forward xLut transform grid points. */ +static void +lutfwdgam_func( + void *pp, /* lutgamctx structure */ + double *out, /* output' value at clut grid point (ie. PCS' value) */ + double *in /* input' value at clut grid point (ie. device' value) */ +) { + lutgamctx *p = (lutgamctx *)pp; + double pcso[3]; /* PCS output value */ + + /* Figure if we are over the ink limit. */ + if ( (p->x->ink.tlimit >= 0.0 || p->x->ink.klimit >= 0.0) + && icxLimitD(p->x, in) > 0.0) { + int i; + double sf; + + /* We are, so use the bracket search to discover a scale */ + /* for the clut input' value that will put us on the ink limit. */ + + for (i = 0; i < p->x->inputChan; i++) + p->in[i] = in[i]; + + if (zbrent(&sf, 0.0, 1.0, 1e-4, icxLimitFind, pp) != 0) { + return; /* Give up */ + } + + /* Compute ink limit value */ + for (i = 0; i < p->x->inputChan; i++) + p->in[i] = sf * in[i]; + + /* Compute the clut output for this clut input */ + p->x->clut(p->x, pcso, p->in); + p->x->output(p->x, pcso, pcso); + p->x->out_abs(p->x, pcso, pcso); + } else { /* No ink limiting */ + /* Convert the clut PCS' values to PCS output values */ + p->x->output(p->x, pcso, out); + p->x->out_abs(p->x, pcso, pcso); + } + + /* Expand the gamut surface with this point */ + p->g->expand(p->g, pcso); + + /* Leave out[] unchanged */ +} + + +/* Function to pass to rspl to create gamut boundary from */ +/* backwards Lut transform. This is called for every node in the */ +/* B2A grid. */ +static void +lutbwdgam_func( + void *pp, /* lutgamctx structure */ + double *out, /* output value */ + double *in /* input value */ +) { + lutgamctx *p = (lutgamctx *)pp; + double devo[MAX_CHAN]; /* Device output value */ + double pcso[3]; /* PCS output value */ + + /* Convert the clut values to device output values */ + p->x->output(p->x, devo, out); /* (Device never uses out_abs()) */ + + /* Convert from device values to PCS values */ + p->flu->lookup(p->flu, pcso, devo); + + /* Expand the gamut surface with this point */ + p->g->expand(p->g, pcso); + + /* Leave out[] unchanged */ +} + +/* Given an xicc lookup object, return a gamut object. */ +/* Note that the PCS must be Lab or Jab */ +/* An icxLuLut type must be icmFwd or icmBwd, */ +/* and for icmFwd, the ink limit (if supplied) will be applied. */ +/* Return NULL on error, check errc+err for reason */ +static gamut *icxLuLutGamut( +icxLuBase *plu, /* this */ +double detail /* gamut detail level, 0.0 = def */ +) { + xicc *p = plu->pp; /* parent xicc */ + icxLuLut *luluto = (icxLuLut *)plu; /* Lookup xLut type object */ + icColorSpaceSignature ins, pcs, outs; + icmLookupFunc func; + icRenderingIntent intent; + double white[3], black[3], kblack[3]; + int inn, outn; + gamut *gam; + + /* get some details */ + plu->spaces(plu, &ins, &inn, &outs, &outn, NULL, &intent, &func, &pcs); + + if (func != icmFwd && func != icmBwd) { + p->errc = 1; + sprintf(p->err,"Creating Gamut surface for anything other than Device <-> PCS is not supported."); + return NULL; + } + + if (pcs != icSigLabData && pcs != icxSigJabData) { + p->errc = 1; + sprintf(p->err,"Creating Gamut surface PCS of other than Lab or Jab is not supported."); + return NULL; + } + + if (func == icmFwd) { + lutgamctx cx; + + cx.g = gam = new_gamut(detail, pcs == icxSigJabData, 0); + cx.x = luluto; + + /* Scan through grid. */ + /* (Note this can give problems for a strange input space - ie. Lab */ + /* and a low grid resolution - ie. 2) */ + luluto->clutTable->scan_rspl( + luluto->clutTable, /* this */ + RSPL_NOFLAGS, /* Combination of flags */ + (void *)&cx, /* Opaque function context */ + lutfwdgam_func /* Function to set from */ + ); + + /* Make sure the white and point goes in too, if it isn't in the grid */ + plu->efv_wh_bk_points(plu, white, NULL, NULL); + gam->expand(gam, white); + + if (detail == 0.0) + detail = 10.0; + + /* If the gamut is more than cursary, add some more detail surface points */ + if (detail < 20.0 || luluto->clutTable->g.mres < 4) { + int res; + DCOUNT(co, MAX_CHAN, inn, 0, 0, 2); + + res = (int)(500.0/detail); /* Establish an appropriate sampling density */ + if (res < 10) + res = 10; + + /* Itterate over all the faces in the device space */ + DC_INIT(co); + while(!DC_DONE(co)) { /* Count through the corners of hyper cube */ + int e, m1, m2; + double in[MAX_CHAN]; + double out[3]; + + for (e = 0; e < inn; e++) { /* Base value */ + in[e] = (double)co[e]; /* Base value */ + in[e] = in[e] * (luluto->inmax[e] - luluto->inmin[e]) + + luluto->inmin[e]; + } + + /* Figure if we are over the ink limit. */ + if ((luluto->ink.tlimit >= 0.0 || luluto->ink.klimit >= 0.0) + && icxLimit(luluto, in) > 0.0) { + DC_INC(co); + continue; /* Skip points over limit */ + } + + /* Scan only device surface */ + for (m1 = 0; m1 < inn; m1++) { /* Choose first coord to scan */ + if (co[m1] != 0) + continue; /* Not at lower corner */ + for (m2 = m1 + 1; m2 < inn; m2++) { /* Choose second coord to scan */ + int x, y; + + if (co[m2] != 0) + continue; /* Not at lower corner */ + + for (x = 0; x < res; x++) { /* step over surface */ + in[m1] = x/(res - 1.0); + in[m1] = in[m1] * (luluto->inmax[m1] - luluto->inmin[m1]) + + luluto->inmin[m1]; + for (y = 0; y < res; y++) { + in[m2] = y/(res - 1.0); + in[m2] = in[m2] * (luluto->inmax[m2] - luluto->inmin[m2]) + + luluto->inmin[m2]; + + /* Figure if we are over the ink limit. */ + if ( (luluto->ink.tlimit >= 0.0 || luluto->ink.klimit >= 0.0) + && icxLimit(luluto, in) > 0.0) { + continue; /* Skip points over limit */ + } + + luluto->lookup((icxLuBase *)luluto, out, in); + gam->expand(gam, out); + } + } + } + } + /* Increment index within block */ + DC_INC(co); + } + } + + /* Now set the cusp points by itterating through colorant 0 & 100% combinations */ + /* If we know what sort of space it is: */ + if (ins == icSigRgbData || ins == icSigCmyData || ins == icSigCmykData) { + DCOUNT(co, 3, 3, 0, 0, 2); + + gam->setcusps(gam, 0, NULL); + DC_INIT(co); + while(!DC_DONE(co)) { + int e; + double in[MAX_CHAN]; + 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]; + in[e] = 0; /* K */ + + /* Always use the device->PCS conversion */ + if (luluto->lookup((icxLuBase *)luluto, out, in) > 1) + error ("%d, %s",p->errc,p->err); + gam->setcusps(gam, 3, out); + } + + DC_INC(co); + } + gam->setcusps(gam, 2, NULL); + + /* Do all ink combinations and hope we can sort it out */ + /* (This may not be smart, since bodgy cusp values will cause gamut mapping to fail...) */ + } else if (ins != icSigXYZData + && ins != icSigLabData + && ins != icSigLuvData + && ins != icSigYxyData) { + DCOUNT(co, MAX_CHAN, inn, 0, 0, 2); + + gam->setcusps(gam, 0, NULL); + DC_INIT(co); + while(!DC_DONE(co)) { + int e; + double in[MAX_CHAN]; + double out[3]; + + for (e = 0; e < inn; e++) { + in[e] = (double)co[e]; + in[e] = in[e] * (luluto->inmax[e] - luluto->inmin[e]) + + luluto->inmin[e]; + } + + /* Figure if we are over the ink limit. */ + if ((luluto->ink.tlimit >= 0.0 || luluto->ink.klimit >= 0.0) + && icxLimit(luluto, in) > 0.0) { + DC_INC(co); + continue; /* Skip points over limit */ + } + + luluto->lookup((icxLuBase *)luluto, out, in); + gam->setcusps(gam, 1, out); + + DC_INC(co); + } + gam->setcusps(gam, 2, NULL); + } + + } else { /* Must be icmBwd */ + lutgamctx cx; + + /* Get an appropriate device to PCS conversion for the fwd conversion */ + /* we use after bwd conversion in lutbwdgam_func() */ + switch (intent) { + /* If it is relative */ + case icmDefaultIntent: /* Shouldn't happen */ + case icPerceptual: + case icRelativeColorimetric: + case icSaturation: + intent = icRelativeColorimetric; /* Choose relative */ + break; + /* If it is absolute */ + case icAbsoluteColorimetric: + case icxAppearance: + case icxAbsAppearance: + break; /* Leave unchanged */ + default: + break; + } + if ((cx.flu = p->get_luobj(p, ICX_CLIP_NEAREST , icmFwd, intent, pcs, icmLuOrdNorm, + &plu->vc, NULL)) == NULL) { + return NULL; /* oops */ + } + + cx.g = gam = new_gamut(detail, pcs == icxSigJabData, 0); + + cx.x = luluto; + + luluto->clutTable->scan_rspl( + luluto->clutTable, /* this */ + RSPL_NOFLAGS, /* Combination of flags */ + (void *)&cx, /* Opaque function context */ + lutbwdgam_func /* Function to set from */ + ); + + /* Now set the cusp points by using the fwd conversion and */ + /* itterating through colorant 0 & 100% combinations. */ + /* If we know what sort of space it is: */ + if (outs == icSigRgbData || outs == icSigCmyData || outs == icSigCmykData) { + DCOUNT(co, 3, 3, 0, 0, 2); + + gam->setcusps(gam, 0, NULL); + DC_INIT(co); + while(!DC_DONE(co)) { + int e; + double in[MAX_CHAN]; + 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]; + in[e] = 0; /* K */ + + /* Always use the device->PCS conversion */ + if (cx.flu->lookup((icxLuBase *)cx.flu, out, in) > 1) + error ("%d, %s",p->errc,p->err); + gam->setcusps(gam, 3, out); + } + + DC_INC(co); + } + gam->setcusps(gam, 2, NULL); + + /* Do all ink combinations and hope we can sort it out */ + /* (This may not be smart, since bodgy cusp values will cause gamut mapping to fail...) */ + } else if (ins != icSigXYZData + && ins != icSigLabData + && ins != icSigLuvData + && ins != icSigYxyData) { + DCOUNT(co, MAX_CHAN, inn, 0, 0, 2); + + gam->setcusps(gam, 0, NULL); + DC_INIT(co); + while(!DC_DONE(co)) { + int e; + double in[MAX_CHAN]; + double out[3]; + + for (e = 0; e < inn; e++) { + in[e] = (double)co[e]; + in[e] = in[e] * (luluto->inmax[e] - luluto->inmin[e]) + + luluto->inmin[e]; + } + + /* Figure if we are over the ink limit. */ + if ((luluto->ink.tlimit >= 0.0 || luluto->ink.klimit >= 0.0) + && icxLimit(luluto, in) > 0.0) { + DC_INC(co); + continue; /* Skip points over limit */ + } + + cx.flu->lookup((icxLuBase *)cx.flu, out, in); + gam->setcusps(gam, 1, out); + + DC_INC(co); + } + gam->setcusps(gam, 2, NULL); + } + cx.flu->del(cx.flu); /* Done with the fwd conversion */ + } + + /* Set the white and black points */ + plu->efv_wh_bk_points(plu, white, black, kblack); + gam->setwb(gam, white, black, kblack); + +//printf("~1 icxLuLutGamut: set black %f %f %f and kblack %f %f %f\n", black[0], black[1], black[2], kblack[0], kblack[1], kblack[2]); + +#ifdef NEVER /* Not sure if this is a good idea ?? */ + /* Since we know this is a profile, force the space and gamut points to be the same */ + gam->getwb(gam, NULL, NULL, white, black, kblack); /* Get the actual gamut white and black points */ + gam->setwb(gam, white, black, kblack); /* Put it back as colorspace one */ +#endif + + return gam; +} + +/* ----------------------------------------------- */ +#ifdef DEBUG +#undef DEBUG /* Limit extent to this file */ +#endif + + + + + + |