diff options
Diffstat (limited to 'xicc/xicc.c')
| -rw-r--r-- | xicc/xicc.c | 3607 |
1 files changed, 3607 insertions, 0 deletions
diff --git a/xicc/xicc.c b/xicc/xicc.c new file mode 100644 index 0000000..9b6d867 --- /dev/null +++ b/xicc/xicc.c @@ -0,0 +1,3607 @@ + +/* + * 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. + * + * Based on the old iccXfm class. + */ + +/* + * This module expands the basic icclib functionality, + * providing more functionality in exercising. + * The implementation for the three different types + * of profile representation, are in their own source files. + */ + +/* + * TTBD: + * Some of the error handling is crude. Shouldn't use + * error(), should return status. + * + */ + +#include <sys/types.h> +#include <string.h> +#include <ctype.h> +#ifdef __sun +#include <unistd.h> +#endif +#if defined(__IBMC__) && defined(_M_IX86) +#include <float.h> +#endif +#include "numlib.h" +#include "counters.h" +#include "plot.h" +#include "../h/sort.h" +#include "xicc.h" /* definitions for this library */ + +#define USE_CAM /* Use CIECAM02 for clipping and gamut mapping, else use Lab */ + +#undef DEBUG /* Plot 1d Luts */ + +#ifdef DEBUG +#include "plot.h" +#endif + +#define MAX_INVSOLN 4 + +static void xicc_del(xicc *p); +icxLuBase * xicc_get_luobj(xicc *p, int flags, icmLookupFunc func, icRenderingIntent intent, + icColorSpaceSignature pcsor, icmLookupOrder order, + icxViewCond *vc, icxInk *ink); +static icxLuBase *xicc_set_luobj(xicc *p, icmLookupFunc func, icRenderingIntent intent, + icmLookupOrder order, int flags, int no, int nobw, cow *points, + icxMatrixModel *skm, + double dispLuminance, double wpscale, double smooth, double avgdev, + icxViewCond *vc, icxInk *ink, xcal *cal, int quality); +static void icxLutSpaces(icxLuBase *p, icColorSpaceSignature *ins, int *inn, + icColorSpaceSignature *outs, int *outn, + icColorSpaceSignature *pcs); +static void icxLuSpaces(icxLuBase *p, icColorSpaceSignature *ins, int *inn, + icColorSpaceSignature *outs, int *outn, + icmLuAlgType *alg, icRenderingIntent *intt, + icmLookupFunc *fnc, icColorSpaceSignature *pcs); +static void icxLu_get_native_ranges (icxLuBase *p, + double *inmin, double *inmax, double *outmin, double *outmax); +static void icxLu_get_ranges (icxLuBase *p, + double *inmin, double *inmax, double *outmin, double *outmax); +static void icxLuEfv_wh_bk_points(icxLuBase *p, double *wht, double *blk, double *kblk); +int xicc_get_viewcond(xicc *p, icxViewCond *vc); + +/* The different profile types are in their own source filesm */ +/* and are included to keep their functions private. (static) */ +#include "xmono.c" +#include "xmatrix.c" +#include "xlut.c" /* New xfit3 in & out optimising based profiles */ +//#include "xlut1.c" /* Old xfit1 device curve based profiles */ + +#ifdef NT /* You'd think there might be some standards.... */ +# ifndef __BORLANDC__ +# define stricmp _stricmp +# endif +#else +# define stricmp strcasecmp +#endif +/* Utilities */ + +/* Return a string description of the given enumeration value */ +const char *icx2str(icmEnumType etype, int enumval) { + + if (etype == icmColorSpaceSignature) { + if (((icColorSpaceSignature)enumval) == icxSigJabData) + return "Jab"; + else if (((icColorSpaceSignature)enumval) == icxSigJChData) + return "JCh"; + else if (((icColorSpaceSignature)enumval) == icxSigLChData) + return "LCh"; + } else if (etype == icmRenderingIntent) { + if (((icRenderingIntent)enumval) == icxAppearance) + return "icxAppearance"; + else if (((icRenderingIntent)enumval) == icxAbsAppearance) + return "icxAbsAppearance"; + else if (((icRenderingIntent)enumval) == icxPerceptualAppearance) + return "icxPerceptualAppearance"; + else if (((icRenderingIntent)enumval) == icxAbsPerceptualAppearance) + return "icxAbsPerceptualAppearance"; + else if (((icRenderingIntent)enumval) == icxSaturationAppearance) + return "icxSaturationAppearance"; + else if (((icRenderingIntent)enumval) == icxAbsSaturationAppearance) + return "icxAbsSaturationAppearance"; + } + return icm2str(etype, enumval); +} + +/* Common xicc stuff */ + +/* Return information about the native lut in/out colorspaces. */ +/* Any pointer may be NULL if value is not to be returned */ +static void +icxLutSpaces( + icxLuBase *p, /* This */ + icColorSpaceSignature *ins, /* Return input color space */ + int *inn, /* Return number of input components */ + icColorSpaceSignature *outs, /* Return output color space */ + int *outn, /* Return number of output components */ + icColorSpaceSignature *pcs /* Return PCS color space */ +) { + p->plu->lutspaces(p->plu, ins, inn, outs, outn, pcs); +} + +/* Return information about the overall lookup in/out colorspaces, */ +/* including allowance for any PCS override. */ +/* Any pointer may be NULL if value is not to be returned */ +static void +icxLuSpaces( + icxLuBase *p, /* This */ + icColorSpaceSignature *ins, /* Return input color space */ + int *inn, /* Return number of input components */ + icColorSpaceSignature *outs, /* Return output color space */ + int *outn, /* Return number of output components */ + icmLuAlgType *alg, /* Return type of lookup algorithm used */ + icRenderingIntent *intt, /* Return the intent implemented */ + icmLookupFunc *fnc, /* Return the profile function being implemented */ + icColorSpaceSignature *pcs /* Return the effective PCS */ +) { + icmLookupFunc function; + icColorSpaceSignature npcs; /* Native PCS */ + + p->plu->spaces(p->plu, NULL, inn, NULL, outn, alg, NULL, &function, &npcs, NULL); + + if (intt != NULL) + *intt = p->intent; + + if (fnc != NULL) + *fnc = function; + + if (ins != NULL) + *ins = p->ins; + + if (outs != NULL) + *outs = p->outs; + + if (pcs != NULL) + *pcs = p->pcs; +} + +/* Return the native (internaly visible) colorspace value ranges */ +static void +icxLu_get_native_ranges ( +icxLuBase *p, +double *inmin, double *inmax, /* Return maximum range of inspace values */ +double *outmin, double *outmax /* Return maximum range of outspace values */ +) { + int i; + if (inmin != NULL) { + for (i = 0; i < p->inputChan; i++) + inmin[i] = p->ninmin[i]; + } + if (inmax != NULL) { + for (i = 0; i < p->inputChan; i++) + inmax[i] = p->ninmax[i]; + } + if (outmin != NULL) { + for (i = 0; i < p->outputChan; i++) + outmin[i] = p->noutmin[i]; + } + if (outmax != NULL) { + for (i = 0; i < p->outputChan; i++) + outmax[i] = p->noutmax[i]; + } +} + +/* Return the effective (externaly visible) colorspace value ranges */ +static void +icxLu_get_ranges ( +icxLuBase *p, +double *inmin, double *inmax, /* Return maximum range of inspace values */ +double *outmin, double *outmax /* Return maximum range of outspace values */ +) { + int i; + if (inmin != NULL) { + for (i = 0; i < p->inputChan; i++) + inmin[i] = p->inmin[i]; + } + if (inmax != NULL) { + for (i = 0; i < p->inputChan; i++) + inmax[i] = p->inmax[i]; + } + if (outmin != NULL) { + for (i = 0; i < p->outputChan; i++) + outmin[i] = p->outmin[i]; + } + if (outmax != NULL) { + for (i = 0; i < p->outputChan; i++) + outmax[i] = p->outmax[i]; + } +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Routine to figure out a suitable black point for CMYK */ + +/* Structure to hold optimisation information */ +typedef struct { + icmLuBase *p; + int kch; /* K channel, -1 if none */ + double tlimit, klimit; /* Ink limit values */ + int inn; /* Number of input channels */ + icColorSpaceSignature outs; /* Output space */ + double p1[3]; /* white pivot point in abs Lab */ + double p2[3]; /* Point on vector towards black in abs Lab */ + double toll; /* Tollerance of black direction */ +} bpfind; + +/* 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 bpfindfunc(void *adata, double pv[]) { + bpfind *b = (bpfind *)adata; + double rv = 0.0; + double Lab[3]; + double lr, ta, tb, terr; /* L ratio, target a, target b, target error */ + double ovr = 0.0; + int e; + + /* Compute amount outside total limit */ + if (b->tlimit >= 0.0) { + double sum; + for (sum = 0.0, e = 0; e < b->inn; e++) + sum += pv[e]; + if (sum > b->tlimit) { + ovr = sum - b->tlimit; +#ifdef DEBUG + printf("~1 total ink ovr = %f\n",ovr); +#endif + } + } + + /* Compute amount outside black limit */ + if (b->klimit >= 0.0 && b->kch >= 0) { + double kval = pv[b->kch] - b->klimit; + if (kval > ovr) { + ovr = kval; +#ifdef DEBUG + printf("~1 black ink ovr = %f\n",ovr); +#endif + } + } + /* Compute amount outside device value limits 0.0 - 1.0 */ + { + double dval; + for (dval = -1.0, e = 0; e < b->inn; e++) { + if (pv[e] < 0.0) { + if (-pv[e] > dval) + dval = -pv[e]; + } else if (pv[e] > 1.0) { + if ((pv[e] - 1.0) > dval) + dval = pv[e] - 1.0; + } + } + if (dval > ovr) + ovr = dval; + } + + /* Compute the Lab value: */ + b->p->lookup(b->p, Lab, pv); + if (b->outs == icSigXYZData) + icmXYZ2Lab(&icmD50, Lab, 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; +} + +/* Try and compute a real black point in XYZ given an iccLu, */ +/* and also return the K only black or the normal black if the device doesn't have K */ +/* black[] will be unchanged if black cannot be computed. */ +/* Note that the black point will be in the space of the Lu */ +/* converted to XYZ, so will have the Lu's intent etc. */ +/* (Note that this is duplicated in xlut.c set_icxLuLut() !!!) */ +static void icxLu_comp_bk_point( +icxLuBase *x, +int gblk, /* If nz, compute black if possible. */ +double *white, /* XYZ Input, used for computing black */ +double *black, /* XYZ Input & Output. Set if gblk NZ and can be computed */ +double *kblack /* XYZ Output. Looked up if possible or set to black[] otherwise */ +) { + icmLuBase *p = x->plu; + icmLuBase *op = p; /* Original icmLu, in case we replace p */ + icc *icco = p->icp; + icmHeader *h = icco->header; + icColorSpaceSignature ins, outs; + int inn, outn; + icmLuAlgType alg; + icRenderingIntent intt; + icmLookupFunc fnc; + icmLookupOrder ord; + int kch = -1; + double dblack[MAX_CHAN]; /* device black value */ + int e; + +#ifdef DEBUG + printf("~1 icxLu_comp_bk_point() called, gblk %d, white = %s, black = %s\n",gblk,icmPdv(3, white),icmPdv(3,black)); +#endif + /* Default return incoming black as K only black */ + kblack[0] = black[0]; + kblack[1] = black[1]; + kblack[2] = black[2]; + + /* Get the effective characteristics of the Lu */ + p->spaces(p, &ins, &inn, &outs, &outn, &alg, &intt, &fnc, NULL, &ord); + + if (fnc == icmBwd) { /* Hmm. We've got PCS to device, and we want device to PCS. */ + + /* Strictly speaking this is a dubious approach, since for a cLut profile */ + /* the B2A table could make the effective white and black points */ + /* anything it likes, and they don't have to match what the corresponding */ + /* A2B table does. In our usage it's probably OK, since we tend */ + /* to use colorimetric B2A */ +#ifdef DEBUG + printf("~1 getting icmFwd\n"); +#endif + if ((p = icco->get_luobj(icco, icmFwd, intt, ins, ord)) == NULL) + error("icxLu_comp_bk_point: assert: getting Fwd Lookup failed!"); + + p->spaces(p, &ins, &inn, &outs, &outn, &alg, &intt, &fnc, NULL, &ord); + } + + if (outs != icSigXYZData && outs != icSigLabData) { + error("icxLu_comp_bk_point: assert: icc Lu output is not XYZ or Lab!, outs = 0x%x, "); + } + +#ifdef DEBUG + printf("~1 icxLu_comp_bk_point called for inn = %d, ins = %s\n", inn, icx2str(icmColorSpaceSignature,ins)); +#endif + + switch (ins) { + + case icSigXYZData: + case icSigLabData: + case icSigLuvData: + case icSigYxyData: +#ifdef DEBUG + printf("~1 Assuming CIE colorspace black is 0.0\n"); +#endif + if (gblk) { + for (e = 0; e < inn; e++) + black[0] = 0.0; + } + kblack[0] = black[0]; + kblack[1] = black[1]; + kblack[2] = black[2]; + return; + + case icSigRgbData: +#ifdef DEBUG + printf("~1 RGB:\n"); +#endif + for (e = 0; e < inn; e++) + dblack[e] = 0.0; + break; + + case icSigGrayData: { /* Could be additive or subtractive */ + double dval[1]; + double minv[3], maxv[3]; +#ifdef DEBUG + printf("~1 Gray:\n"); +#endif + /* Check out 0 and 100% colorant */ + dval[0] = 0.0; + p->lookup(p, minv, dval); + if (outs == icSigXYZData) + icmXYZ2Lab(&icmD50, minv, minv); + dval[0] = 1.0; + p->lookup(p, maxv, dval); + if (outs == icSigXYZData) + icmXYZ2Lab(&icmD50, maxv, maxv); + + if (minv[0] < maxv[0]) + dblack[0] = 0.0; + else + dblack[0] = 1.0; + } + break; + + case icSigCmyData: + for (e = 0; e < inn; e++) + dblack[e] = 1.0; + break; + + case icSigCmykData: +#ifdef DEBUG + printf("~1 CMYK:\n"); +#endif + kch = 3; + dblack[0] = 0.0; + dblack[1] = 0.0; + dblack[2] = 0.0; + dblack[3] = 1.0; + if (alg == icmLutType) { + icxLuLut *pp = (icxLuLut *)x; + + if (pp->ink.tlimit >= 0.0) + dblack[kch] = pp->ink.tlimit; + }; + break; + + /* Use a heursistic. */ + /* This duplicates code in icxGetLimits() :-( */ + /* Colorant guessing should go in icclib ? */ + case icSig2colorData: + case icSig3colorData: + case icSig4colorData: + case icSig5colorData: + case icSig6colorData: + case icSig7colorData: + case icSig8colorData: + case icSig9colorData: + case icSig10colorData: + case icSig11colorData: + case icSig12colorData: + case icSig13colorData: + case icSig14colorData: + case icSig15colorData: + case icSigMch5Data: + case icSigMch6Data: + case icSigMch7Data: + case icSigMch8Data: { + double dval[MAX_CHAN]; + double ncval[3]; + double cvals[MAX_CHAN][3]; + int nlighter, ndarker; + + /* Decide if the colorspace is additive or subtractive */ +#ifdef DEBUG + printf("~1 N channel:\n"); +#endif + + /* First the no colorant value */ + for (e = 0; e < inn; e++) + dval[e] = 0.0; + p->lookup(p, ncval, dval); + if (outs == icSigXYZData) + icmXYZ2Lab(&icmD50, ncval, ncval); + + /* Then all the colorants */ + nlighter = ndarker = 0; + for (e = 0; e < inn; e++) { + dval[e] = 1.0; + p->lookup(p, cvals[e], dval); + if (outs == icSigXYZData) + icmXYZ2Lab(&icmD50, cvals[e], cvals[e]); + dval[e] = 0.0; + if (fabs(cvals[e][0] - ncval[0]) > 5.0) { + if (cvals[e][0] > ncval[0]) + nlighter++; + else + ndarker++; + } + } + if (ndarker == 0 && nlighter > 0) { /* Assume additive */ + for (e = 0; e < inn; e++) + dblack[e] = 0.0; +#ifdef DEBUG + printf("~1 N channel is additive:\n"); +#endif + + } else if (ndarker > 0 && nlighter == 0) { /* Assume subtractive. */ + double pbk[3] = { 0.0,0.0,0.0 }; /* Perfect black */ + double smd = 1e10; /* Smallest distance */ + +#ifdef DEBUG + printf("~1 N channel is subtractive:\n"); +#endif + /* See if we can guess the black channel */ + for (e = 0; e < inn; e++) { + double tt; + tt = icmNorm33sq(pbk, cvals[e]); + if (tt < smd) { + smd = tt; + kch = e; + } + } + /* See if the black seems sane */ + if (cvals[kch][0] > 40.0 + || fabs(cvals[kch][1]) > 10.0 + || fabs(cvals[kch][2]) > 10.0) { + if (p != op) + p->del(p); +#ifdef DEBUG + printf("~1 black doesn't look sanem so assume nothing\n"); +#endif + return; /* Assume nothing */ + } + + /* Chosen kch as black */ + for (e = 0; e < inn; e++) + dblack[e] = 0.0; + dblack[kch] = 1.0; + if (alg == icmLutType) { + icxLuLut *pp = (icxLuLut *)x; + + if (pp->ink.tlimit >= 0.0) + dblack[kch] = pp->ink.tlimit; + }; +#ifdef DEBUG + printf("~1 N channel K = chan %d\n",kch); +#endif + } else { + if (p != op) + p->del(p); +#ifdef DEBUG + printf("~1 can't figure if additive or subtractive, so assume nothing\n"); +#endif + return; /* Assume nothing */ + } + } + break; + + default: +#ifdef DEBUG + printf("~1 unhandled colorspace, so assume nothing\n"); +#endif + if (p != op) + p->del(p); + return; /* Don't do anything */ + } + + /* Lookup the K only value */ + if (kch >= 0) { + p->lookup(p, kblack, dblack); + + /* We always return XYZ */ + if (outs == icSigLabData) + icmLab2XYZ(&icmD50, kblack, kblack); + } + + if (gblk == 0) { /* That's all we have to do */ +#ifdef DEBUG + printf("~1 gblk == 0, so only return kblack\n"); +#endif + if (p != op) + p->del(p); + return; + } + + /* Lookup the device black or K only value as a default */ + p->lookup(p, black, dblack); /* May be XYZ or Lab */ +#ifdef DEBUG + printf("~1 Got default lu black %f %f %f, kch = %d\n", black[0],black[1],black[2],kch); +#endif + + /* !!! 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 ?? */ + + if (kch >= 0) { /* The space is subtractive with a K channel. */ + /* If XICC_NEUTRAL_CMYK_BLACK then locate the darkest */ + /* CMYK within limits with the same chromaticity as the white point, */ + /* otherwise locate the device value within the ink limits that is */ + /* in the direction of the K channel */ + bpfind bfs; /* Callback context */ + double sr[MXDO]; /* search radius */ + double tt[MXDO]; /* Temporary */ + double rs0[MXDO], rs1[MXDO]; /* Random start candidates */ + int trial; + double brv; + + /* Setup callback function context */ + bfs.p = p; + bfs.inn = inn; + bfs.outs = outs; + + bfs.kch = kch; + bfs.tlimit = -1.0; + bfs.klimit = -1.0; + bfs.toll = XICC_BLACK_POINT_TOLL; + + if (alg == icmLutType) { + icxLuLut *pp = (icxLuLut *)x; + + pp->kch = kch; + bfs.tlimit = pp->ink.tlimit; + bfs.klimit = pp->ink.klimit; +#ifdef DEBUG + printf("~1 tlimit = %f, klimit = %f\n",bfs.tlimit,bfs.klimit); +#endif + }; + +#ifdef XICC_NEUTRAL_CMYK_BLACK +#ifdef DEBUG + printf("~1 Searching for neutral black\n"); +#endif + /* white has been given to us in XYZ */ + icmXYZ2Lab(&icmD50, bfs.p1, white); /* pivot Lab */ + icmCpy3(bfs.p2, white); /* temp white XYZ */ + icmScale3(bfs.p2, bfs.p2, 0.02); /* Scale white XYZ towards 0,0,0 */ + icmXYZ2Lab(&icmD50, bfs.p2, bfs.p2); /* Convert black direction to Lab */ + +#else /* Use K directin black */ +#ifdef DEBUG + printf("~1 Searching for K direction black\n"); +#endif + icmXYZ2Lab(&icmD50, bfs.p1, white); /* Pivot */ + + /* Now figure abs Lab value of K only, as the direction */ + /* to use for the rich black. */ + for (e = 0; e < inn; e++) + dblack[e] = 0.0; + if (bfs.klimit < 0.0) + dblack[kch] = 1.0; + else + dblack[kch] = bfs.klimit; /* K value */ + + p->lookup(p, black, dblack); + + if (outs == icSigXYZData) { + icmXYZ2Lab(&icmD50, bfs.p2, black); /* K direction */ + } else { + icmAry2Ary(bfs.p2, black); + } +#endif + +#ifdef DEBUG + printf("~1 Lab pivot %f %f %f, Lab K direction %f %f %f\n",bfs.p1[0],bfs.p1[1],bfs.p1[2],bfs.p2[0],bfs.p2[1],bfs.p2[2]); +#endif + /* Start with the K only as the current best value */ + brv = bpfindfunc((void *)&bfs, dblack); +#ifdef DEBUG + printf("~1 initial brv for K only = %f\n",brv); +#endif + + /* Set the random start 0 location as 000K */ + /* and the random start 1 location as CMY0 */ + { + double tt; + + for (e = 0; e < inn; e++) + dblack[e] = rs0[e] = 0.0; + if (bfs.klimit < 0.0) + dblack[kch] = rs0[kch] = 1.0; + else + dblack[kch] = rs0[kch] = bfs.klimit; /* K value */ + + if (bfs.tlimit < 0.0) + tt = 1.0; + else + tt = bfs.tlimit/(inn - 1.0); + for (e = 0; e < inn; e++) + rs1[e] = tt; + 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 (e = 0; e < inn; e++) { + tt[e] = rs0[e]; + sr[e] = 0.1; + } + + } else { + /* Base is random between 000K and CMY0: */ + if (trial < 100) { + rv = d_rand(0.0, 1.0); + for (e = 0; e < inn; e++) { + tt[e] = rv * rs0[e] + (1.0 - rv) * rs1[e]; + sr[e] = 0.1; + } + /* Base on current best */ + } else { + for (e = 0; e < inn; e++) { + tt[e] = dblack[e]; + sr[e] = 0.1; + } + } + + /* Then add random start offset */ + for (rv = 0.0, e = 0; e < inn; e++) { + tt[e] += d_rand(-0.5, 0.5); + if (tt[e] < 0.0) + tt[e] = 0.0; + else if (tt[e] > 1.0) + tt[e] = 1.0; + } + } + + /* Clip black */ + if (bfs.klimit >= 0.0 && tt[kch] > bfs.klimit) + tt[kch] = bfs.klimit; + + /* Compute amount outside total limit */ + if (bfs.tlimit >= 0.0) { + for (rv = 0.0, e = 0; e < inn; e++) + rv += tt[e]; + + if (rv > bfs.tlimit) { + rv /= (double)inn; + for (e = 0; e < inn; e++) + tt[e] -= rv; + } + } + + if (powell(&rv, inn, tt, sr, 0.000001, 1000, bpfindfunc, + (void *)&bfs, NULL, NULL) == 0) { +#ifdef DEBUG + printf("~1 trial %d, rv %f bp %f %f %f %f\n",trial,rv,tt[0],tt[1],tt[2],tt[3]); +#endif + if (rv < brv) { +#ifdef DEBUG + printf("~1 new best\n"); +#endif + brv = rv; + for (e = 0; e < inn; e++) + dblack[e] = tt[e]; + } + } + } + if (brv > 1000.0) + error("icxLu_comp_bk_point: Black point powell failed"); + + for (e = 0; e < inn; e++) { /* Make sure device values are in range */ + if (dblack[e] < 0.0) + dblack[e] = 0.0; + else if (dblack[e] > 1.0) + dblack[e] = 1.0; + } + /* Now have device black in dblack[] */ +#ifdef DEBUG + printf("~1 got device black %f %f %f %f\n",dblack[0], dblack[1], dblack[2], dblack[3]); +#endif + + p->lookup(p, black, dblack); /* Convert to PCS */ + } + + if (p != op) + p->del(p); + + /* We always return XYZ */ + if (outs == icSigLabData) + icmLab2XYZ(&icmD50, black, black); + +#ifdef DEBUG + printf("~1 returning %f %f %f\n", black[0], black[1], black[2]); +#endif + + return; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - */ + +/* Return the media white and black points */ +/* in the xlu effective PCS colorspace. Pointers may be NULL. */ +/* (ie. these will be relative values for relative intent etc.) */ +static void icxLuEfv_wh_bk_points( +icxLuBase *p, +double *wht, +double *blk, +double *kblk /* K only black */ +) { + double white[3], black[3], kblack[3]; + + /* Get the Lu PCS converted to XYZ icc black and white points in XYZ */ + if (p->plu->lu_wh_bk_points(p->plu, white, black)) { + /* Black point is assumed. We should determine one instead. */ + /* Lookup K only black too */ + icxLu_comp_bk_point(p, 1, white, black, kblack); + + } else { + /* Lookup a possible K only black */ + icxLu_comp_bk_point(p, 0, white, black, kblack); + } + +//printf("~1 white %f %f %f, black %f %f %f, kblack %f %f %f\n",white[0],white[1],white[2],black[0],black[1],black[2],kblack[0],kblack[1],kblack[2]); + /* Convert to possibl xicc override PCS */ + switch (p->pcs) { + case icSigXYZData: + break; /* Don't have to do anyting */ + case icSigLabData: + icmXYZ2Lab(&icmD50, white, white); /* Convert from XYZ to Lab */ + icmXYZ2Lab(&icmD50, black, black); + icmXYZ2Lab(&icmD50, kblack, kblack); + break; + case icxSigJabData: + p->cam->XYZ_to_cam(p->cam, white, white); /* Convert from XYZ to Jab */ + p->cam->XYZ_to_cam(p->cam, black, black); + p->cam->XYZ_to_cam(p->cam, kblack, kblack); + break; + default: + break; + } + +//printf("~1 icxLuEfv_wh_bk_points: pcsor %s White %f %f %f, Black %f %f %f\n", icx2str(icmColorSpaceSignature,p->pcs), white[0], white[1], white[2], black[0], black[1], black[2]); + if (wht != NULL) { + wht[0] = white[0]; + wht[1] = white[1]; + wht[2] = white[2]; + } + + if (blk != NULL) { + blk[0] = black[0]; + blk[1] = black[1]; + blk[2] = black[2]; + } + + if (kblk != NULL) { + kblk[0] = kblack[0]; + kblk[1] = kblack[1]; + kblk[2] = kblack[2]; + } +} + +/* Create an instance of an xicc object */ +xicc *new_xicc( +icc *picc /* icc we are expanding */ +) { + xicc *p; + if ((p = (xicc *) calloc(1,sizeof(xicc))) == NULL) + return NULL; + p->pp = picc; + p->del = xicc_del; + p->get_luobj = xicc_get_luobj; + p->set_luobj = xicc_set_luobj; + p->get_viewcond = xicc_get_viewcond; + + /* Create an xcal if there is the right tag in the profile */ + p->cal = xiccReadCalTag(p->pp); + p->nodel_cal = 0; /* We created it, we will delete it */ + + return p; +} + +/* Do away with the xicc (but not the icc!) */ +static void xicc_del( +xicc *p +) { + if (p->cal != NULL && p->nodel_cal == 0) + p->cal->del(p->cal); + free (p); +} + + +/* Return an expanded lookup object, initialised */ +/* from the icc. */ +/* Return NULL on error, check errc+err for reason. */ +/* Set the pcsor & intent to consistent and values if */ +/* Jab and/or icxAppearance has been requested. */ +/* Create 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 this affects the individual */ +/* conversion steps, which will all talk the native PCS (unless merged). */ +icxLuBase *xicc_get_luobj( +xicc *p, /* this */ +int flags, /* clip, merge flags */ +icmLookupFunc func, /* Functionality */ +icRenderingIntent intent, /* Intent */ +icColorSpaceSignature pcsor,/* PCS override (0 = def) */ +icmLookupOrder order, /* Search Order */ +icxViewCond *vc, /* Viewing Condition (may be NULL if pcsor is not CIECAM) */ +icxInk *ink /* inking details (NULL for default) */ +) { + icmLuBase *plu; + icxLuBase *xplu; + icmLuAlgType alg; + icRenderingIntent n_intent = intent; /* Native Intent to request */ + icColorSpaceSignature n_pcs = icmSigDefaultData; /* Native PCS to request */ + +//printf("~1 xicc_get_luobj got intent %s and pcsor %s\n",icx2str(icmRenderingIntent,intent),icx2str(icmColorSpaceSignature,pcsor)); + + /* Ensure that appropriate PCS is slected for an appearance intent */ + if (intent == icxAppearance + || intent == icxAbsAppearance + || intent == icxPerceptualAppearance + || intent == icxAbsPerceptualAppearance + || intent == icxSaturationAppearance + || intent == icxAbsSaturationAppearance) { + pcsor = icxSigJabData; + + /* Translate non-Jab intents to the equivalent appearance "intent" if pcsor == Jab. */ + /* This is how we get these when the UI's don't list all the apperances intents, */ + /* we select the analogous non-apperance intent with pcsor = Jab. */ + /* Note that Abs/non-abs selects between Apperance and AbsAppearance. */ + } else if (pcsor == icxSigJabData) { + if (intent == icRelativeColorimetric) + intent = icxAppearance; + else if (intent == icAbsoluteColorimetric) + intent = icxAbsAppearance; + else if (intent == icPerceptual) + intent = icxPerceptualAppearance; + else if (intent == icmAbsolutePerceptual) + intent = icxAbsPerceptualAppearance; + else if (intent == icSaturation) + intent = icxSaturationAppearance; + else if (intent == icmAbsoluteSaturation) + intent = icxAbsSaturationAppearance; + else + intent = icxAppearance; + } + + /* Translate intent asked for into intent needed in icclib */ + if (intent == icxAppearance + || intent == icxAbsAppearance) + n_intent = icAbsoluteColorimetric; + else if (intent == icxPerceptualAppearance + || intent == icxAbsPerceptualAppearance) + n_intent = icmAbsolutePerceptual; + else if (intent == icxSaturationAppearance + || intent == icxAbsSaturationAppearance) + n_intent = icmAbsoluteSaturation; + + if (pcsor != icmSigDefaultData) + n_pcs = pcsor; /* There is an icclib override */ + + if (pcsor == icxSigJabData) /* xicc override */ + n_pcs = icSigXYZData; /* Translate to XYZ */ + +//printf("~1 xicc_get_luobj processed intent %s and pcsor %s\n",icx2str(icmRenderingIntent,intent),icx2str(icmColorSpaceSignature,pcsor)); +//printf("~1 xicc_get_luobj icclib intent %s and pcsor %s\n",icx2str(icmRenderingIntent,n_intent),icx2str(icmColorSpaceSignature,n_pcs)); + /* Get icclib lookup object */ + if ((plu = p->pp->get_luobj(p->pp, func, n_intent, n_pcs, order)) == NULL) { + p->errc = p->pp->errc; /* Copy error */ + strcpy(p->err, p->pp->err); + return NULL; + } + + /* Figure out what the algorithm is */ + plu->spaces(plu, NULL, NULL, NULL, NULL, &alg, NULL, NULL, &n_pcs, NULL); + + /* make sure its "Abs CAM" */ + if (vc!= NULL + && (intent == icxAbsAppearance + || intent == icxAbsPerceptualAppearance + || intent == icxAbsSaturationAppearance)) { /* make sure its "Abs CAM" */ + /* Set white point and flare color to D50 */ + /* (Hmm. This doesn't match what happens within collink with absolute intent!!) */ + vc->Wxyz[0] = icmD50.X/icmD50.Y; + vc->Wxyz[1] = icmD50.Y/icmD50.Y; // Normalise white reference to Y = 1 ? + vc->Wxyz[2] = icmD50.Z/icmD50.Y; + + vc->Fxyz[0] = icmD50.X; + vc->Fxyz[1] = icmD50.Y; + vc->Fxyz[2] = icmD50.Z; + } + + /* Call xiccLu wrapper creation */ + switch (alg) { + case icmMonoFwdType: + xplu = new_icxLuMono(p, flags, plu, func, intent, pcsor, vc, 0); + break; + case icmMonoBwdType: + xplu = new_icxLuMono(p, flags, plu, func, intent, pcsor, vc, 1); + break; + case icmMatrixFwdType: + xplu = new_icxLuMatrix(p, flags, plu, func, intent, pcsor, vc, 0); + break; + case icmMatrixBwdType: + xplu = new_icxLuMatrix(p, flags, plu, func, intent, pcsor, vc, 1); + break; + case icmLutType: + xplu = new_icxLuLut(p, flags, plu, func, intent, pcsor, vc, ink); + break; + default: + xplu = NULL; + break; + } + + return xplu; +} + + +/* Return an expanded lookup object, initialised */ +/* from the icc, and then overwritten by a conversion */ +/* created from the supplied scattered data points. */ +/* The Lut is assumed to be a device -> native PCS profile. */ +/* If the SET_WHITE and/or SET_BLACK flags are set, */ +/* discover the white/black point, set it in the icc, */ +/* and make the Lut relative to them. */ +/* Return NULL on error, check errc+err for reason */ +static icxLuBase *xicc_set_luobj( +xicc *p, /* this */ +icmLookupFunc func, /* Functionality */ +icRenderingIntent intent, /* Intent */ +icmLookupOrder order, /* Search Order */ +int flags, /* white/black point, verbose flags etc. */ +int no, /* Number of points */ +int nobw, /* Number of points to look for white & black patches in */ +cow *points, /* Array of input points in target PCS space */ +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 input white point is to be scaled */ +double smooth, /* RSPL smoothing factor, -ve if raw */ +double avgdev, /* reading Average Deviation as a proportion of the input range */ +icxViewCond *vc, /* Viewing Condition (NULL if not using CAM) */ +icxInk *ink, /* inking details (NULL for default) */ +xcal *cal, /* Optional cal, will override any existing (not deleted with xicc)*/ +int quality /* Quality metric, 0..3 */ +) { + icmLuBase *plu; + icxLuBase *xplu = NULL; + icmLuAlgType alg; + + if (cal != NULL) { + if (p->cal != NULL && p->nodel_cal == 0) + p->cal->del(p->cal); + p->cal = cal; + p->nodel_cal = 1; /* We were given it, so don't delete it */ + } + + if (func != icmFwd) { + p->errc = 1; + sprintf(p->err,"Can only create Device->PCS profiles from scattered data."); + xplu = NULL; + return xplu; + } + + /* Get icclib lookup object */ + if ((plu = p->pp->get_luobj(p->pp, func, intent, 0, order)) == NULL) { + p->errc = p->pp->errc; /* Copy error */ + strcpy(p->err, p->pp->err); + return NULL; + } + + /* Figure out what the algorithm is */ + plu->spaces(plu, NULL, NULL, NULL, NULL, &alg, NULL, NULL, NULL, NULL); + + /* Call xiccLu wrapper creation */ + switch (alg) { + case icmMonoFwdType: + p->errc = 1; + sprintf(p->err,"Setting Monochrome Fwd profile from scattered data not supported."); + plu->del(plu); + xplu = NULL; /* Not supported yet */ + break; + + case icmMatrixFwdType: + if (smooth < 0.0) + smooth = -smooth; + xplu = set_icxLuMatrix(p, plu, flags, no, nobw, points, skm, dispLuminance, wpscale, quality, smooth); + break; + + case icmLutType: + /* ~~~ Should add check that it is a fwd profile ~~~ */ + xplu = set_icxLuLut(p, plu, func, intent, flags, no, nobw, points, skm, dispLuminance, wpscale, smooth, avgdev, vc, ink, quality); + break; + + default: + break; + } + + return xplu; +} + +/* ------------------------------------------------------ */ +/* Viewing Condition Parameter stuff */ + +#ifdef NEVER /* Not currently used */ + +/* Guess viewing parameters from the technology signature */ +static void guess_from_techsig( +icTechnologySignature tsig, +double *Ybp +) { +double Yb = -1.0; + + switch (tsig) { + /* These are all inputing either a representation of */ + /* a natural scene captured on another medium, or are assuming */ + /* that the medium is the original. A _good_ system would */ + /* let the user indicate which is the case. */ + case icSigReflectiveScanner: + case icSigFilmScanner: + Yb = 0.2; + break; + + /* Direct scene to value devices. */ + case icSigDigitalCamera: + case icSigVideoCamera: + Yb = 0.2; + break; + + /* Emmisive displays. */ + /* We could try tweaking the white point on the assumption */ + /* that the viewer will be adapted to a combination of both */ + /* the CRT white point, and the ambient light. */ + case icSigVideoMonitor: + case icSigCRTDisplay: + case icSigPMDisplay: + case icSigAMDisplay: + Yb = 0.2; + break; + + /* Photo CD has its own viewing definitions */ + /* (It represents original scene colors) */ + case icSigPhotoCD: + Yb = 0.2; + break; + + /* Projection devices, either direct, or */ + /* via another intermediate medium. */ + case icSigProjectionTelevision: + Yb = 0.1; /* Assume darkened room, little background */ + break; + case icSigFilmWriter: + Yb = 0.0; /* Assume a dark room - no background */ + break; + + /* Printed media devices. */ + case icSigInkJetPrinter: + case icSigThermalWaxPrinter: + case icSigElectrophotographicPrinter: + case icSigElectrostaticPrinter: + case icSigDyeSublimationPrinter: + case icSigPhotographicPaperPrinter: + case icSigPhotoImageSetter: + case icSigGravure: + case icSigOffsetLithography: + case icSigSilkscreen: + case icSigFlexography: + Yb = 0.2; + break; + + default: + Yb = 0.2; + } + + if (Ybp != NULL) + *Ybp = Yb; +} + +#endif /* NEVER */ + + +/* See if we can read or guess the viewing conditions */ +/* for an ICC profile. */ +/* Return value 0 if it is well defined */ +/* Return value 1 if it is a guess */ +/* Return value 2 if it is not possible/appropriate */ +int xicc_get_viewcond( +xicc *p, /* Expanded profile we're working with */ +icxViewCond *vc /* Viewing parameters to return */ +) { + icc *pp = p->pp; /* Base ICC */ + + /* Numbers we're trying to find */ + ViewingCondition Ev = vc_none; + double Wxyz[3] = {-1.0, -1.0, -1.0}; /* Adapting white color */ + double La = -1.0; /* Adapting luminance */ + double Ixyz[3] = {-1.0, -1.0, -1.0}; /* Illuminant color */ + double Li = -1.0; /* Illuminant luminance */ + double Lb = -1.0; /* Backgrount luminance */ + double Yb = -1.0; /* Background relative luminance to Lv */ + double Lve = -1.0; /* Emissive device image luminance */ + double Lvr = -1.0; /* Reflective device image luminance */ + double Lv = -1.0; /* device image luminance */ + double Yf = -1.0; /* Flare relative luminance to Lv */ + double Fxyz[3] = {-1.0, -1.0, -1.0}; /* Flare color */ + icTechnologySignature tsig = icMaxEnumTechnology; /* Technology Signature */ + icProfileClassSignature devc = icMaxEnumClass; + int trans = -1; /* Set to 0 if not transparency, 1 if it is */ + + /* Collect all the information we can find */ + + /* Emmisive devices image white luminance */ + { + icmXYZArray *luminanceTag; + + if ((luminanceTag = (icmXYZArray *)pp->read_tag(pp, icSigLuminanceTag)) != NULL + && luminanceTag->ttype == icSigXYZType && luminanceTag->size >= 1) { + Lve = luminanceTag->data[0].Y; /* Copy structure */ + } + } + + /* Flare: */ + { + icmMeasurement *ro; + + if ((ro = (icmMeasurement *)pp->read_tag(pp, icSigMeasurementTag)) != NULL + && ro->ttype == icSigMeasurementType) { + + Yf = ro->flare; + /* ro->illuminant ie D50, D65, D93, A etc. */ + } + } + + /* Media White Point */ + { + icmXYZArray *whitePointTag; + + if ((whitePointTag = (icmXYZArray *)pp->read_tag(pp, icSigMediaWhitePointTag)) != NULL + && whitePointTag->ttype == icSigXYZType && whitePointTag->size >= 1) { + Wxyz[0] = whitePointTag->data[0].X; + Wxyz[1] = whitePointTag->data[0].Y; + Wxyz[2] = whitePointTag->data[0].Z; + } + } + + /* ViewingConditions: */ + { + icmViewingConditions *ro; + + if ((ro = (icmViewingConditions *)pp->read_tag(pp, icSigViewingConditionsTag)) != NULL + && ro->ttype == icSigViewingConditionsType) { + + /* ro->illuminant.X */ + /* ro->illuminant.Z */ + + Li = ro->illuminant.Y; + + /* Reflect illuminant off the media white */ + Lvr = Li * Wxyz[1]; + + /* Illuminant color */ + Ixyz[0] = ro->illuminant.X/ro->illuminant.Y; + Ixyz[1] = 1.0; + Ixyz[2] = ro->illuminant.Z/ro->illuminant.Y; + + /* Assume ICC surround is CICAM97 background */ + /* ro->surround.X */ + /* ro->surround.Z */ + La = ro->surround.Y; + + /* ro->stdIlluminant ie D50, D65, D93, A etc. */ + } + } + + /* Stuff we might need */ + + /* Technology: */ + { + icmSignature *ro; + + /* Try and read the tag from the file */ + if ((ro = (icmSignature *)pp->read_tag(pp, icSigTechnologyTag)) != NULL + && ro->ttype != icSigSignatureType) { + + tsig = ro->sig; + } + } + + devc = pp->header->deviceClass; /* Type of profile */ + if (devc == icSigLinkClass + || devc == icSigAbstractClass + || devc == icSigColorSpaceClass + || devc == icSigNamedColorClass) + return 2; + + /* + icSigInputClass + icSigDisplayClass + icSigOutputClass + */ + + if ((pp->header->flags & icTransparency) != 0) + trans = 1; + else + trans = 0; + + + /* figure Lv if we have the information */ + if (Lve >= 0.0) + Lv = Lve; /* Emmisive image white luminance */ + else + Lv = Lvr; /* Reflectance image white luminance */ + + /* Fudge the technology signature */ + if (tsig == icMaxEnumTechnology) { + if (devc == icSigDisplayClass) + tsig = icSigCRTDisplay; + } + +#ifndef NEVER + printf("Enumeration = %d\n", Ev); + printf("Viewing Conditions:\n"); + printf("White adaptation color %f %f %f\n",Wxyz[0], Wxyz[1], Wxyz[2]); + printf("Adapting Luminance La = %f\n",La); + printf("Illuminant color %f %f %f\n",Ixyz[0], Ixyz[1], Ixyz[2]); + printf("Illuminant Luminance Li = %f\n",Li); + printf("Background Luminance Lb = %f\n",Lb); + printf("Relative Background Yb = %f\n",Yb); + printf("Emissive Image White Lve = %f\n",Lve); + printf("Reflective Image White Lvr = %f\n",Lvr); + printf("Device Image White Lv = %f\n",Lv); + printf("Relative Flare Yf = %f\n",Yf); + printf("Flare color %f %f %f\n",Fxyz[0], Fxyz[1], Fxyz[2]); + printf("Technology = %s\n",tag2str(tsig)); + printf("deviceClass = %s\n",tag2str(devc)); + printf("Transparency = %d\n",trans); +#endif + + /* See if the viewing conditions are completely defined as ICC can do it */ + if (Wxyz[0] >= 0.0 && Wxyz[1] >= 0.0 && Wxyz[2] >= 0.0 + && La >= 0.0 + && Yb >= 0.0 + && Lv >= 0.0 + && Yf >= 0.0 + && Fxyz[0] >= 0.0 && Fxyz[1] >= 0.0 && Fxyz[2] >= 0.0) { + + vc->Ev = vc_none; + vc->Wxyz[0] = Wxyz[0]; + vc->Wxyz[1] = Wxyz[1]; + vc->Wxyz[2] = Wxyz[2]; + vc->La = La; + vc->Yb = Yb; + vc->Lv = Lv; + vc->Yf = Yf; + vc->Fxyz[0] = Fxyz[0]; + vc->Fxyz[1] = Fxyz[1]; + vc->Fxyz[2] = Fxyz[2]; + return 0; + } + + /* Hmm. We didn't get all the info an ICC can contain. */ + /* We will try to guess some reasonable defaults */ + + /* Have we at least got an adaptation white point ? */ + if (Wxyz[0] < 0.0 || Wxyz[1] < 0.0 || Wxyz[2] < 0.0) + return 2; /* No */ + + /* Have we got the technology ? */ + if (tsig == icMaxEnumTechnology) + return 2; /* Hopeless */ + + /* Guess from the technology */ + switch (tsig) { + + /* This is inputing either a representation of */ + /* a natural scene captured on another a print medium, or */ + /* are is assuming that the medium is the original. */ + /* We will assume that the print is the original. */ + case icSigReflectiveScanner: + { + if (La < 0.0) /* No adapting luminance */ + La = 34.0; /* Use a practical print evaluation number */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.2; /* Assume grey world */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_average; /* Assume average viewing conditions */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.01; /* Assume 1% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + + /* This is inputing either a representation of */ + /* a natural scene captured on another a photo medium, or */ + /* are is assuming that the medium is the original. */ + /* We will assume a compromise media original, natural scene */ + case icSigFilmScanner: + { + if (La < 0.0) /* No adapting luminance */ + La = 50.0; /* Use bright indoors, dull outdoors */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.2; /* Assume grey world */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_average; /* Assume average viewing conditions */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.005; /* Assume 0.5% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + + /* Direct scene to value devices. */ + case icSigDigitalCamera: + case icSigVideoCamera: + { + if (La < 0.0) /* No adapting luminance */ + La = 110.0; /* Use very bright indoors, usual outdoors */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.2; /* Assume grey world */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_average; /* Assume average viewing conditions */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.0; /* Assume 0% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + + /* Emmisive displays. */ + /* Assume a video monitor is in a darker environment than a CRT */ + case icSigVideoMonitor: + { + if (La < 0.0) /* No adapting luminance */ + La = 4.0; /* Darkened work environment */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.2; /* Assume grey world */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_dim; /* Assume dim viewing conditions */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.01; /* Assume 1% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + + + /* Assume a typical work environment */ + case icSigCRTDisplay: + case icSigPMDisplay: + case icSigAMDisplay: + { + if (La < 0.0) /* No adapting luminance */ + La = 33.0; /* Typical work environment */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.2; /* Assume grey world */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_average; /* Assume average viewing conditions */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.02; /* Assume 2% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + + /* Photo CD has its own viewing definitions */ + /* (It represents original scene colors) */ + case icSigPhotoCD: + { + if (La < 0.0) /* No adapting luminance */ + La = 320.0; /* Bright outdoors */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.2; /* Assume grey world */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_average; /* Assume average viewing conditions */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.00; /* Assume 0% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + + /* Projection devices, either direct, or */ + /* via another intermediate medium. */ + /* Assume darkened room, little background */ + case icSigProjectionTelevision: + { + if (La < 0.0) /* No adapting luminance */ + La = 7.0; /* Dark environment */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.1; /* Assume little background */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_dim; /* Dim environment */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.01; /* Assume 1% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + /* Assume very darkened room, no background */ + case icSigFilmWriter: + { + if (La < 0.0) /* No adapting luminance */ + La = 7.0; /* Dark environment */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.0; /* Assume no background */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_dark; /* Dark environment */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.01; /* Assume 1% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + + /* Printed media devices. */ + /* Assume a normal print viewing environment */ + case icSigInkJetPrinter: + case icSigThermalWaxPrinter: + case icSigElectrophotographicPrinter: + case icSigElectrostaticPrinter: + case icSigDyeSublimationPrinter: + case icSigPhotographicPaperPrinter: + case icSigPhotoImageSetter: + case icSigGravure: + case icSigOffsetLithography: + case icSigSilkscreen: + case icSigFlexography: + { + if (La < 0.0) /* No adapting luminance */ + La = 40.0; /* Use a practical print evaluation number */ + if (Yb < 0.0) /* No background relative luminance */ + Yb = 0.2; /* Assume grey world */ + if (Lv < 0.0) /* No device image luminance */ + Ev = vc_average; /* Assume average viewing conditions */ + if (Yf < 0.0) /* No flare figure */ + Yf = 0.01; /* Assume 1% flare */ + if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0) /* No flare color */ + Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2]; + break; + } + + default: + { + return 2; + } + } + + return 1; +} + +/* Write our viewing conditions to the underlying ICC profile, */ +/* using a private tag. */ +void xicc_set_viewcond( +xicc *p, /* Expanded profile we're working with */ +icxViewCond *vc /* Viewing parameters to return */ +) { + //icc *pp = p->pp; /* Base ICC */ + + // ~~1 Not implemented yet +} + + + +/* Return an enumerated viewing condition */ +/* Return enumeration if OK, -999 if there is no such enumeration. */ +/* xicc may be NULL if just the description is wanted, */ +/* or an explicit white point is provided. */ +int xicc_enum_viewcond( +xicc *p, /* Expanded profile to get white point (May be NULL if desc NZ) */ +icxViewCond *vc, /* Viewing parameters to return, May be NULL if desc is nz */ +int no, /* Enumeration to return, -1 for default, -2 for none */ +char *as, /* String alias to number, NULL if none */ +int desc, /* NZ - Just return a description of this enumeration in vc */ +double *wp /* Provide white point if xicc is NULL */ +) { + + if (desc == 0) { /* We're setting the viewing condition */ + icc *pp; /* Base ICC */ + icmXYZArray *whitePointTag; + + if (vc == NULL) + return -999; + + if (p == NULL) { + if (wp == NULL) + return -999; + vc->Wxyz[0] = wp[0]; + vc->Wxyz[1] = wp[1]; + vc->Wxyz[2] = wp[2]; + } else { + + pp = p->pp; + if ((whitePointTag = (icmXYZArray *)pp->read_tag(pp, icSigMediaWhitePointTag)) != NULL + && whitePointTag->ttype == icSigXYZType && whitePointTag->size >= 1) { + vc->Wxyz[0] = whitePointTag->data[0].X; + vc->Wxyz[1] = whitePointTag->data[0].Y; + vc->Wxyz[2] = whitePointTag->data[0].Z; + } else { + if (wp == NULL) { + sprintf(p->err,"Enum VC: Failed to read Media White point"); + p->errc = 2; + return -999; + } + vc->Wxyz[0] = wp[0]; + vc->Wxyz[1] = wp[1]; + vc->Wxyz[2] = wp[2]; + } + } + + /* Set a default flare color */ + vc->Fxyz[0] = vc->Wxyz[0]; + vc->Fxyz[1] = vc->Wxyz[1]; + vc->Fxyz[2] = vc->Wxyz[2]; + } + + /* + + Typical adapting field luminances and white luminance in reflective setup: + + E = illuminance in Lux + Lv = White luminance assuming 100% reflectance + La = Adapting field luminance in cd/m^2, assuming 20% reflectance from surround + + E La Lv Condition + 11 0.7 4 Twilight + 32 2 10 Subdued indoor lighting + 64 4 20 Less than typical office light; sometimes recommended for + display-only workplaces (sRGB) + 350 22 111 Typical Office (sRGB annex D) + 500 32 160 Practical print evaluationa (ISO-3664 P2) + 1000 64 318 Good Print evaluation (CIE 116-1995) + 1000 64 318 Television Studio lighting + 1000 64 318 Overcast Outdoors + 2000 127 637 Critical print evaluation (ISO-3664 P1) + 10000 637 3183 Typical outdoors, full daylight + 50000 3185 15915 Bright summers day + + */ + + if (no == -1 + || (as != NULL && stricmp(as,"d") == 0)) { + + no = -1; + if (vc != NULL) { + vc->desc = " d - Default Viewing Condition"; + vc->Ev = vc_average; /* Average viewing conditions */ + vc->La = 50.0; /* Practical to Good lighting */ + vc->Lv = 250.0; /* Average viewing conditions ratio */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.01; /* 1% flare */ + } + } + else if (no == 0 + || (as != NULL && stricmp(as,"pp") == 0)) { + + no = 0; + if (vc != NULL) { + vc->desc = " pp - Practical Reflection Print (ISO-3664 P2)"; + vc->Ev = vc_average; /* Average viewing conditions */ + vc->La = 32.0; /* Use a practical print evaluation number */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.01; /* 1% flare */ + } + } + else if (no == 1 + || (as != NULL && stricmp(as,"pe") == 0)) { + + no = 1; + if (vc != NULL) { + vc->desc = " pe - Print evaluation environment (CIE 116-1995)"; + vc->Ev = vc_average; /* Average viewing conditions */ + vc->La = 64.0; /* Good */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.01; /* 1% flare */ + } + } + else if (no == 2 + || (as != NULL && stricmp(as,"pc") == 0)) { + + no = 2; + if (vc != NULL) { + vc->desc = " pc - Critical print evaluation environment (ISO-3664 P1)"; + vc->Ev = vc_average; /* Average viewing conditions */ + vc->La = 127.0; /* Critical */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.01; /* 1% flare */ + } + } + else if (no == 3 + || (as != NULL && stricmp(as,"mt") == 0)) { + + no = 3; + if (vc != NULL) { + vc->desc = " mt - Monitor in typical work environment"; + vc->Ev = vc_average; /* Average viewing conditions */ + vc->La = 22.0; /* Typical work environment */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.02; /* 2% flare */ + } + } + else if (no == 4 + || (as != NULL && stricmp(as,"mb") == 0)) { + + no = 4; + if (vc != NULL) { + vc->desc = " mb - Bright monitor in bright work environment"; + vc->Ev = vc_average; /* Average viewing conditions */ + vc->La = 42.0; /* Bright work environment */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.02; /* 2% flare */ + } + } + else if (no == 5 + || (as != NULL && stricmp(as,"md") == 0)) { + + no = 5; + if (vc != NULL) { + vc->desc = " md - Monitor in darkened work environment"; + vc->Ev = vc_dim; /* Dim viewing conditions */ + vc->La = 4.0; /* Darkened work environment */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.01; /* 1% flare */ + } + } + else if (no == 6 + || (as != NULL && stricmp(as,"jm") == 0)) { + + no = 6; + if (vc != NULL) { + vc->desc = " jm - Projector in dim environment"; + vc->Ev = vc_dim; /* Dim viewing conditions */ + vc->La = 10.0; /* Adaptation is from display */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.01; /* 1% flare */ + } + } + else if (no == 7 + || (as != NULL && stricmp(as,"jd") == 0)) { + + no = 7; + if (vc != NULL) { + vc->desc = " jd - Projector in dark environment"; + vc->Ev = vc_dark; /* Dark viewing conditions */ + vc->La = 10.0; /* Adaptation is from display */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.01; /* 1% flare ? */ + } + } + else if (no == 8 + || (as != NULL && stricmp(as,"pcd") == 0)) { + + no = 8; + if (vc != NULL) { + vc->desc = "pcd - Photo CD - original scene outdoors"; + vc->Ev = vc_average; /* Average viewing conditions */ + vc->La = 320.0; /* Typical outdoors, 1600 cd/m^2 */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.00; /* 0% flare */ + } + } + else if (no == 9 + || (as != NULL && stricmp(as,"ob") == 0)) { + + no = 9; + if (vc != NULL) { + vc->desc = " ob - Original scene - Bright Outdoors"; + vc->Ev = vc_average; /* Average viewing conditions */ + vc->La = 2000.0; /* Bright Outdoors */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.00; /* 0% flare */ + } + } + else if (no == 10 + || (as != NULL && stricmp(as,"cx") == 0)) { + + no = 10; + if (vc != NULL) { + vc->desc = " cx - Cut Sheet Transparencies on a viewing box"; + vc->Ev = vc_cut_sheet; /* Cut sheet viewing conditions */ + vc->La = 53.0; /* Dim, adapted to slide ? */ + vc->Yb = 0.2; /* Grey world */ + vc->Yf = 0.01; /* 1% flare ? */ + } + } + else { + if (p != NULL) { + sprintf(p->err,"Enum VC: Unrecognised enumeration %d",no); + p->errc = 1; + } + return -999; + } + + return no; +} + +/* Debug: dump a Viewing Condition to standard out */ +void xicc_dump_viewcond( +icxViewCond *vc +) { + printf("Viewing Condition:\n"); + if (vc->Ev == vc_dark) + printf(" Surround to Image: Dark\n"); + else if (vc->Ev == vc_dim) + printf(" Surround to Image: Dim\n"); + else if (vc->Ev == vc_average) + printf(" Surround to Image: Average\n"); + else if (vc->Ev == vc_cut_sheet) + printf(" Transparency on Light box\n"); + printf(" Adapted white = %f %f %f\n",vc->Wxyz[0], vc->Wxyz[1], vc->Wxyz[2]); + printf(" Adapted luminance = %f cd/m^2\n",vc->La); + printf(" Background to image ratio = %f\n",vc->Yb); + if (vc->Ev == vc_none) + printf(" Image luminance = %f cd/m^2\n",vc->Lv); + printf(" Flare to image ratio = %f\n",vc->Yf); + printf(" Flare color = %f %f %f\n",vc->Fxyz[0], vc->Fxyz[1], vc->Fxyz[2]); +} + + +/* Debug: dump an Inking setup to standard out */ +void xicc_dump_inking(icxInk *ik) { + printf("Inking settings:\n"); + if (ik->tlimit < 0.0) + printf("No total limit\n"); + else + printf("Total limit = %f%%\n",ik->tlimit * 100.0); + + if (ik->klimit < 0.0) + printf("No black limit\n"); + else + printf("Black limit = %f%%\n",ik->klimit * 100.0); + + if (ik->KonlyLmin) + printf("K only black as locus Lmin\n"); + else + printf("Normal black as locus Lmin\n"); + + if (ik->k_rule == icxKvalue) { + printf("Inking rule is a fixed K target\n"); + } if (ik->k_rule == icxKlocus) { + printf("Inking rule is a fixed locus target\n"); + } if (ik->k_rule == icxKluma5 || ik->k_rule == icxKluma5k) { + if (ik->k_rule == icxKluma5) + printf("Inking rule is a 5 parameter locus function of L\n"); + else + printf("Inking rule is a 5 parameter K function of L\n"); + printf("Ksmth = %f\n",ik->c.Ksmth); + printf("Kskew = %f\n",ik->c.Kskew); + printf("Kstle = %f\n",ik->c.Kstle); + printf("Kstpo = %f\n",ik->c.Kstpo); + printf("Kenpo = %f\n",ik->c.Kenpo); + printf("Kenle = %f\n",ik->c.Kenle); + printf("Kshap = %f\n",ik->c.Kshap); + } if (ik->k_rule == icxKl5l || ik->k_rule == icxKl5lk) { + if (ik->k_rule == icxKl5l) + printf("Inking rule is a 2x5 parameter locus function of L and K aux\n"); + else + printf("Inking rule is a 2x5 parameter K function of L and K aux\n"); + printf("Min Ksmth = %f\n",ik->c.Ksmth); + printf("Min Kskew = %f\n",ik->c.Kskew); + printf("Min Kstle = %f\n",ik->c.Kstle); + printf("Min Kstpo = %f\n",ik->c.Kstpo); + printf("Min Kenpo = %f\n",ik->c.Kenpo); + printf("Min Kenle = %f\n",ik->c.Kenle); + printf("Min Kshap = %f\n",ik->c.Kshap); + printf("Max Ksmth = %f\n",ik->x.Ksmth); + printf("Max Kskew = %f\n",ik->x.Kskew); + printf("Max Kstle = %f\n",ik->x.Kstle); + printf("Max Kstpo = %f\n",ik->x.Kstpo); + printf("Max Kenpo = %f\n",ik->x.Kenpo); + printf("Max Kenle = %f\n",ik->x.Kenle); + printf("Max Kshap = %f\n",ik->x.Kshap); + } +} + +/* ------------------------------------------------------ */ +/* Gamut Mapping Intent stuff */ + +/* Return an enumerated gamut mapping intent */ +/* Return enumeration if OK, icxIllegalGMIntent if there is no such enumeration. */ +int xicc_enum_gmapintent( +icxGMappingIntent *gmi, /* Gamut Mapping parameters to return */ +int no, /* Enumeration selected, icxNoGMIntent for none */ +char *as /* Alias string selector, NULL for none */ +) { +#ifdef USE_CAM + int colccas = 0x2; /* Use cas clipping for colorimetric style intents */ + int perccas = 0x1; /* Use cas for perceptual style intents */ +#else + int colccas = 0x0; /* Use Lab for colorimetric style intents */ + int perccas = 0x0; /* Use Lab for perceptual style intents */ + fprintf(stderr,"!!!!!! Warning, USE_CAM is off in xicc.c !!!!!!\n"); +#endif + + /* Assert default if no guidance given */ + if (no == icxNoGMIntent && as == NULL) + no = icxDefaultGMIntent; + + if (no == 0 + || no == icxAbsoluteGMIntent + || (as != NULL && stricmp(as,"a") == 0)) { + /* Map Absolute appearance space Jab to Jab and clip out of gamut */ + no = 0; + gmi->as = "a"; + gmi->desc = " a - Absolute Colorimetric (in Jab) [ICC Absolute Colorimetric]"; + gmi->icci = icAbsoluteColorimetric; + gmi->usecas = colccas; /* Use absolute appearance space */ + gmi->usemap = 0; /* Don't use gamut mapping */ + gmi->greymf = 0.0; + gmi->glumwcpf = 0.0; + gmi->glumwexf = 0.0; + gmi->glumbcpf = 0.0; + gmi->glumbexf = 0.0; + gmi->glumknf = 0.0; + gmi->gamcpf = 0.0; + gmi->gamexf = 0.0; + gmi->gamcknf = 0.0; + gmi->gamxknf = 0.0; + gmi->gampwf = 0.0; + gmi->gamswf = 0.0; + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 1 + || (as != NULL && stricmp(as,"aw") == 0)) { + + /* I'm not sure how often this intent is useful. It's less likely than */ + /* I though that a printer white point won't fit within the gamut */ + /* of a display profile, since the display white always has Y = 1.0, */ + /* and no paper has better than about 95% reflectance. */ + /* Perhaps it may be more useful for targeting printer profiles ? */ + + /* Map Absolute Jab to Jab and scale source to avoid clipping the white point */ + no = 1; + gmi->as = "aw"; + gmi->desc = "aw - Absolute Colorimetric (in Jab) with scaling to fit white point"; + gmi->icci = icAbsoluteColorimetric; + gmi->usecas = 0x100 | colccas; /* Absolute Appearance space with scaling */ + /* to avoid clipping the source white point */ + gmi->usemap = 0; /* Don't use gamut mapping */ + gmi->greymf = 0.0; + gmi->glumwcpf = 0.0; + gmi->glumwexf = 0.0; + gmi->glumbcpf = 0.0; + gmi->glumbexf = 0.0; + gmi->glumknf = 0.0; + gmi->gamcpf = 0.0; + gmi->gamexf = 0.0; + gmi->gamcknf = 0.0; + gmi->gamxknf = 0.0; + gmi->gampwf = 0.0; + gmi->gamswf = 0.0; + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 2 + || (as != NULL && stricmp(as,"aa") == 0)) { + + /* Map appearance space Jab to Jab and clip out of gamut */ + no = 2; + gmi->as = "aa"; + gmi->desc = "aa - Absolute Appearance"; + gmi->icci = icRelativeColorimetric; + gmi->usecas = perccas; /* Appearance space */ + gmi->usemap = 0; /* Don't use gamut mapping */ + gmi->greymf = 0.0; + gmi->glumwcpf = 0.0; + gmi->glumwexf = 0.0; + gmi->glumbcpf = 0.0; + gmi->glumbexf = 0.0; + gmi->glumknf = 0.0; + gmi->gamcpf = 0.0; + gmi->gamexf = 0.0; + gmi->gamcknf = 0.0; + gmi->gamxknf = 0.0; + gmi->gampwf = 0.0; + gmi->gamswf = 0.0; + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 3 + || no == icxRelativeGMIntent + || (as != NULL && stricmp(as,"r") == 0)) { + + /* Align neutral axes and linearly map white point, then */ + /* map appearance space Jab to Jab and clip out of gamut */ + no = 3; + gmi->as = "r"; + gmi->desc = " r - White Point Matched Appearance [ICC Relative Colorimetric]"; + gmi->icci = icRelativeColorimetric; + gmi->usecas = perccas; /* Appearance space */ + gmi->usemap = 1; /* Use gamut mapping */ + gmi->greymf = 1.0; /* Fully align grey axis */ + gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */ + gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */ + gmi->glumbcpf = 0.0; /* No compression at black end */ + gmi->glumbexf = 0.0; /* No expansion at black end */ + gmi->glumknf = 0.0; + gmi->gamcpf = 0.0; + gmi->gamexf = 0.0; + gmi->gamcknf = 0.0; + gmi->gamxknf = 0.0; + gmi->gampwf = 0.0; + gmi->gamswf = 0.0; + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 4 + || (as != NULL && stricmp(as,"la") == 0)) { + + /* Align neutral axes and linearly map white and black points, then */ + /* map appearance space Jab to Jab and clip out of gamut */ + no = 4; + gmi->as = "la"; + gmi->desc = "la - Luminance axis matched Appearance"; + gmi->icci = icRelativeColorimetric; + gmi->usecas = perccas; /* Appearance space */ + gmi->usemap = 1; /* Use gamut mapping */ + gmi->greymf = 1.0; /* Fully align grey axis */ + gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */ + gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */ + gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */ + gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */ + gmi->glumknf = 0.0; /* No knee on grey mapping */ + gmi->gamcpf = 0.0; /* No gamut compression */ + gmi->gamexf = 0.0; /* No gamut expansion */ + gmi->gamcknf = 0.0; /* No knee in gamut compress */ + gmi->gamxknf = 0.0; /* No knee in gamut expand */ + gmi->gampwf = 0.0; /* No Perceptual surface weighting factor */ + gmi->gamswf = 0.0; /* No Saturation surface weighting factor */ + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 5 + || no == icxDefaultGMIntent + || no == icxPerceptualGMIntent + || (as != NULL && stricmp(as,"p") == 0)) { + + /* Align neutral axes and perceptually map white and black points, */ + /* perceptually compress out of gamut and map appearance space Jab to Jab. */ + no = 5; + gmi->as = "p"; + gmi->desc = " p - Perceptual (Preferred) (Default) [ICC Perceptual]"; + gmi->icci = icPerceptual; + gmi->usecas = perccas; /* Appearance space */ + gmi->usemap = 1; /* Use gamut mapping */ + gmi->greymf = 1.0; /* Fully align grey axis */ + gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */ + gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */ + gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */ + gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */ + gmi->glumknf = 1.0; /* Sigma knee in grey compress/expand */ + gmi->gamcpf = 1.0; /* Full gamut compression */ + gmi->gamexf = 0.0; /* No gamut expansion */ + gmi->gamcknf = 0.8; /* High Sigma knee in gamut compress */ + gmi->gamxknf = 0.0; /* No knee in gamut expand */ + gmi->gampwf = 1.0; /* Full Perceptual surface weighting factor */ + gmi->gamswf = 0.0; /* No Saturation surface weighting factor */ + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 6 + || (as != NULL && stricmp(as,"pa") == 0)) { + + /* Don't align neutral axes, but perceptually compress out of gamut */ + /* and map appearance space Jab to Jab. */ + no = 5; + gmi->as = "pa"; + gmi->desc = "pa - Perceptual Apperance "; + gmi->icci = icPerceptual; + gmi->usecas = perccas; /* Appearance space */ + gmi->usemap = 1; /* Use gamut mapping */ + gmi->greymf = 0.0; /* Don't align grey axis */ + gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */ + gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */ + gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */ + gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */ + gmi->glumknf = 1.0; /* Sigma knee in grey compress/expand */ + gmi->gamcpf = 1.0; /* Full gamut compression */ + gmi->gamexf = 0.0; /* No gamut expansion */ + gmi->gamcknf = 0.8; /* High Sigma knee in gamut compress */ + gmi->gamxknf = 0.0; /* No knee in gamut expand */ + gmi->gampwf = 1.0; /* Full Perceptual surface weighting factor */ + gmi->gamswf = 0.0; /* No Saturation surface weighting factor */ + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 7 + || (as != NULL && stricmp(as,"ms") == 0)) { + + /* Align neutral axes and perceptually map white and black points, */ + /* perceptually compress and expand to match gamuts and map Jab to Jab. */ + no = 6; + gmi->as = "ms"; + gmi->desc = "ms - Saturation"; + gmi->icci = icSaturation; + gmi->usecas = perccas; /* Appearance space */ + gmi->usemap = 1; /* Use gamut mapping */ + gmi->greymf = 1.0; /* Fully align grey axis */ + gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */ + gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */ + gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */ + gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */ + gmi->glumknf = 1.0; /* Sigma knee in grey compress/expand */ + gmi->gamcpf = 1.0; /* Full gamut compression */ + gmi->gamexf = 1.0; /* Full gamut expansion */ + gmi->gamcknf = 1.0; /* High Sigma knee in gamut compress/expand */ + gmi->gamxknf = 0.4; /* Moderate Sigma knee in gamut compress/expand */ + gmi->gampwf = 0.2; /* Slight perceptual surface weighting factor */ + gmi->gamswf = 0.8; /* Most saturation surface weighting factor */ + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 8 + || no == icxSaturationGMIntent + || (as != NULL && stricmp(as,"s") == 0)) { + + /* Same as "ms" but enhance saturation */ + no = 7; + gmi->as = "s"; + gmi->desc = " s - Enhanced Saturation [ICC Saturation]"; + gmi->icci = icSaturation; + gmi->usecas = perccas; /* Appearance space */ + gmi->usemap = 1; /* Use gamut mapping */ + gmi->greymf = 1.0; /* Fully align grey axis */ + gmi->glumwcpf = 1.0; /* Fully compress grey axis at white end */ + gmi->glumwexf = 1.0; /* Fully expand grey axis at white end */ + gmi->glumbcpf = 1.0; /* Fully compress grey axis at black end */ + gmi->glumbexf = 1.0; /* Fully expand grey axis at black end */ + gmi->glumknf = 1.0; /* Sigma knee in grey compress/expand */ + gmi->gamcpf = 1.0; /* Full gamut compression */ + gmi->gamexf = 1.0; /* Full gamut expansion */ + gmi->gamcknf = 1.0; /* High sigma knee in gamut compress */ + gmi->gamxknf = 0.5; /* Moderate sigma knee in gamut expand */ + gmi->gampwf = 0.0; /* No Perceptual surface weighting factor */ + gmi->gamswf = 1.0; /* Full Saturation surface weighting factor */ + gmi->satenh = 0.9; /* Medium saturation enhancement */ + } + else if (no == 9 + || (as != NULL && stricmp(as,"al") == 0)) { + + /* Map absolute L*a*b* to L*a*b* and clip out of gamut */ + no = 8; + gmi->as = "al"; + gmi->desc = "al - Absolute Colorimetric (Lab)"; + gmi->icci = icAbsoluteColorimetric; + gmi->usecas = 0x0; /* Don't use appearance space, use L*a*b* */ + gmi->usemap = 0; /* Don't use gamut mapping */ + gmi->greymf = 0.0; + gmi->glumwcpf = 0.0; + gmi->glumwexf = 0.0; + gmi->glumbcpf = 0.0; + gmi->glumbexf = 0.0; + gmi->glumknf = 0.0; + gmi->gamcpf = 0.0; + gmi->gamexf = 0.0; + gmi->gamcknf = 0.0; + gmi->gamxknf = 0.0; + gmi->gampwf = 0.0; + gmi->gamswf = 0.0; + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else if (no == 10 + || (as != NULL && stricmp(as,"rl") == 0)) { + + /* Align neutral axes and linearly map white point, then */ + /* map L*a*b* to L*a*b* and clip out of gamut */ + no = 3; + gmi->as = "rl"; + gmi->desc = "rl - White Point Matched Appearance (Lab)"; + gmi->icci = icRelativeColorimetric; + gmi->usecas = 0x0; /* Don't use appearance space, use L*a*b* */ + gmi->usemap = 1; /* Use gamut mapping */ + gmi->greymf = 1.0; /* And linearly map white point */ + gmi->glumwcpf = 1.0; + gmi->glumwexf = 1.0; + gmi->glumbcpf = 0.0; + gmi->glumbexf = 0.0; + gmi->glumknf = 0.0; + gmi->gamcpf = 0.0; + gmi->gamexf = 0.0; + gmi->gamcknf = 0.0; + gmi->gamxknf = 0.0; + gmi->gampwf = 0.0; + gmi->gamswf = 0.0; + gmi->satenh = 0.0; /* No saturation enhancement */ + } + else { /* icxIllegalGMIntent */ + return icxIllegalGMIntent; + } + + return no; +} + + +/* Debug: dump a Gamut Mapping specification */ +void xicc_dump_gmi( +icxGMappingIntent *gmi /* Gamut Mapping parameters to return */ +) { + printf(" Gamut Mapping Specification:\n"); + if (gmi->desc != NULL) + printf(" Description = '%s'\n",gmi->desc); + printf(" Closest ICC intent = '%s'\n",icm2str(icmRenderingIntent,gmi->icci)); + + if ((gmi->usecas & 0xff) == 0) + printf(" Not using Color Apperance Space\n"); + else if ((gmi->usecas & 0xff) == 1) + printf(" Using Color Apperance Space\n"); + else if ((gmi->usecas & 0xff) == 2) + printf(" Using Absolute Color Apperance Space\n"); + + if ((gmi->usecas & 0x100) != 0) + printf(" Scaling source to avoid white point clipping\n"); + + if (gmi->usemap == 0) + printf(" Not using Mapping\n"); + else { + printf(" Using Mapping with parameters:\n"); + printf(" Grey axis alignment factor %f\n", gmi->greymf); + printf(" Grey axis white compression factor %f\n", gmi->glumwcpf); + printf(" Grey axis white expansion factor %f\n", gmi->glumwexf); + printf(" Grey axis black compression factor %f\n", gmi->glumbcpf); + printf(" Grey axis black expansion factor %f\n", gmi->glumbexf); + printf(" Grey axis knee factor %f\n", gmi->glumknf); + printf(" Gamut compression factor %f\n", gmi->gamcpf); + printf(" Gamut expansion factor %f\n", gmi->gamexf); + printf(" Gamut compression knee factor %f\n", gmi->gamcknf); + printf(" Gamut expansion knee factor %f\n", gmi->gamxknf); + printf(" Gamut Perceptual mapping weighting factor %f\n", gmi->gampwf); + printf(" Gamut Saturation mapping weighting factor %f\n", gmi->gamswf); + printf(" Saturation enhancement factor %f\n", gmi->satenh); + } +} + +/* ------------------------------------------------------ */ +/* Turn xicc xcal into limit calibration callback */ + +/* Given an icc profile, try and create an xcal */ +/* Return NULL on error or no cal */ +xcal *xiccReadCalTag(icc *p) { + xcal *cal = NULL; + icTagSignature sig = icmMakeTag('t','a','r','g'); + icmText *ro; + int oi, tab; + +//printf("~1 about to look for CAL in profile\n"); + if ((ro = (icmText *)p->read_tag(p, sig)) != NULL) { + cgatsFile *cgf; + cgats *icg; + + if (ro->ttype != icSigTextType) + return NULL; + +//printf("~1 found 'targ' tag\n"); + if ((icg = new_cgats()) == NULL) { + return NULL; + } + if ((cgf = new_cgatsFileMem(ro->data, ro->size)) != NULL) { + icg->add_other(icg, "CTI3"); + oi = icg->add_other(icg, "CAL"); + +//printf("~1 created cgats object from 'targ' tag\n"); + if (icg->read(icg, cgf) == 0) { + + for (tab = 0; tab < icg->ntables; tab++) { + if (icg->t[tab].tt == tt_other && icg->t[tab].oi == oi) { + break; + } + } + if (tab < icg->ntables) { +//printf("~1 found CAL table\n"); + + if ((cal = new_xcal()) == NULL) { + icg->del(icg); + cgf->del(cgf); + return NULL; + } + if (cal->read_cgats(cal, icg, tab, "'targ' tag") != 0) { +#ifdef DEBUG + printf("read_cgats on cal tag failed\n"); +#endif + cal->del(cal); + cal = NULL; + } +//else printf("~1 read CAL and creaded xcal object OK\n"); + } + } + cgf->del(cgf); + } + icg->del(icg); + } + return cal; +} + +/* A callback that uses an xcal, that can be used with icc get_tac */ +void xiccCalCallback(void *cntx, double *out, double *in) { + xcal *cal = (xcal *)cntx; + + cal->interp(cal, out, in); +} + +/* ---------------------------------------------- */ + +/* Utility function - given an open icc profile, */ +/* guess which channel is the black. */ +/* Return -1 if there is no black channel or it can't be guessed */ +int icxGuessBlackChan(icc *p) { + int kch = -1; + + switch (p->header->colorSpace) { + case icSigCmykData: + kch = 3; + break; + + /* Use a heuristic to detect the black channel. */ + /* This duplicates code in icxLu_comp_bk_point() :-( */ + /* Colorant guessing should go in icclib ? */ + case icSig2colorData: + case icSig3colorData: + case icSig4colorData: + case icSig5colorData: + case icSig6colorData: + case icSig7colorData: + case icSig8colorData: + case icSig9colorData: + case icSig10colorData: + case icSig11colorData: + case icSig12colorData: + case icSig13colorData: + case icSig14colorData: + case icSig15colorData: + case icSigMch5Data: + case icSigMch6Data: + case icSigMch7Data: + case icSigMch8Data: { + icmLuBase *lu; + double dval[MAX_CHAN]; + double ncval[3]; + double cvals[MAX_CHAN][3]; + int inn, e, nlighter, ndarker; + + /* Grab a lookup object */ + if ((lu = p->get_luobj(p, icmFwd, icRelativeColorimetric, icSigLabData, icmLuOrdNorm)) == NULL) + error("icxGetLimits: assert: getting Fwd Lookup failed!"); + + lu->spaces(lu, NULL, &inn, NULL, NULL, NULL, NULL, NULL, NULL, NULL); + + /* Decide if the colorspace is aditive or subtractive */ + + /* First the no colorant value */ + for (e = 0; e < inn; e++) + dval[e] = 0.0; + lu->lookup(lu, ncval, dval); + + /* Then all the colorants */ + nlighter = ndarker = 0; + for (e = 0; e < inn; e++) { + dval[e] = 1.0; + lu->lookup(lu, cvals[e], dval); + dval[e] = 0.0; + if (fabs(cvals[e][0] - ncval[0]) > 5.0) { + if (cvals[e][0] > ncval[0]) + nlighter++; + else + ndarker++; + } + } + + if (ndarker > 0 && nlighter == 0) { /* Assume subtractive. */ + double pbk[3] = { 0.0,0.0,0.0 }; /* Perfect black */ + double smd = 1e10; /* Smallest distance */ + + /* Guess the black channel */ + for (e = 0; e < inn; e++) { + double tt; + tt = icmNorm33sq(pbk, cvals[e]); + if (tt < smd) { + smd = tt; + kch = e; + } + } + /* See if the black seems sane */ + if (cvals[kch][0] > 40.0 + || fabs(cvals[kch][1]) > 10.0 + || fabs(cvals[kch][2]) > 10.0) { + kch = -1; + } + } + lu->del(lu); + } + break; + + default: + break; + } + + return kch; +} + +/* Utility function - given an open icc profile, */ +/* estmate the total ink limit and black ink limit. */ +/* Note that this is rather rough, because ICC profiles */ +/* don't have a tag for this information, and ICC profiles */ +/* don't have any straightforward way of identifying particular */ +/* color channels for > 4 color. */ +/* If there are no limits, or they are not discoverable or */ +/* applicable, return values of -1.0 */ + +void icxGetLimits( +xicc *xp, +double *tlimit, +double *klimit +) { + icc *p = xp->pp; + int nch; + double max[MAX_CHAN]; /* Max of each channel */ + double total; + + total = p->get_tac(p, max, xp->cal != NULL ? xiccCalCallback : NULL, (void *)xp->cal); + + if (total < 0.0) { /* Not valid */ + if (tlimit != NULL) + *tlimit = -1.0; + if (klimit != NULL) + *klimit = -1.0; + return; + } + + nch = icmCSSig2nchan(p->header->colorSpace); + + /* No effective limit */ + if (tlimit != NULL) { + if (total >= (double)nch) { + *tlimit = -1.0; + } else { + *tlimit = total; + } + } + + if (klimit != NULL) { + int kch; + + kch = icxGuessBlackChan(p); + + if (kch < 0 || max[kch] >= 1.0) { + *klimit = -1.0; + } else { + *klimit = max[kch]; + } + } +} + +/* Replace a non-set limit (ie. < 0.0) with the heuristic from */ +/* the given profile. */ +void icxDefaultLimits( +xicc *xp, +double *tlout, +double tlin, +double *klout, +double klin +) { + if (tlin < 0.0 || klin < 0.0) { + double tl, kl; + + icxGetLimits(xp, &tl, &kl); + + if (tlin < 0.0) + tlin = tl; + + if (klin < 0.0) + klin = kl; + } + + if (tlout != NULL) + *tlout = tlin; + + if (klout != NULL) + *klout = klin; +} + +/* Structure to hold optimisation information */ +typedef struct { + xcal *cal; + double ilimit; + double uilimit; +} ulimctx; + +/* Callback to find equivalent underlying total limit */ +/* and try and maximize it while remaining within gamut */ +static double ulimitfunc(void *cntx, double pv[]) { + ulimctx *cx = (ulimctx *)cntx; + xcal *cal = cx->cal; + int devchan = cal->devchan; + int i; + double dv, odv; + double og = 0.0, rv = 0.0; + + double usum = 0.0, sum = 0.0; + + /* Comute calibrated sum of channels except last */ + for (i = 0; i < (devchan-1); i++) { + double dv = pv[i]; /* Underlying (pre-calibration) device value */ + usum += dv; /* Underlying sum */ + if (dv < 0.0) { + og += -dv; + dv = 0.0; + } else if (dv > 1.0) { + og += dv - 1.0; + dv = 1.0; + } else + dv = cal->interp_ch(cal, i, dv); /* Calibrated device value */ + sum += dv; /* Calibrated device sum */ + } + /* Compute the omitted channel value */ + dv = cx->ilimit - sum; /* Omitted calibrated device value */ + if (dv < 0.0) { + og += -dv; + dv = 0.0; + } else if (dv > 1.0) { + og += dv - 1.0; + dv = 1.0; + } else + dv = cal->inv_interp_ch(cal, i, dv); /* Omitted underlying device value */ + usum += dv; /* Underlying sum */ + cx->uilimit = usum; + + rv = 10000.0 * og - usum; /* Penalize out of gamut, maximize underlying sum */ + +//printf("~1 returning %f from %f %f %f %f\n",rv,pv[0],pv[1],pv[2],dv); + return rv; +} + +/* Given a calibrated total ink limit and an xcal, return the */ +/* equivalent underlying (pre-calibration) total ink limit. */ +/* This is the maximum equivalent, that makes sure that */ +/* the calibrated limit is met or exceeded. */ +double icxMaxUnderlyingLimit(xcal *cal, double ilimit) { + ulimctx cx; + int i; + double dv[MAX_CHAN]; + double sr[MAX_CHAN]; + double rv; /* Residual value */ + + if (cal->devchan <= 1) { + return cal->inv_interp_ch(cal, 0, ilimit); + } + + cx.cal = cal; + cx.ilimit = ilimit; + + for (i = 0; i < (cal->devchan-1); i++) { + sr[i] = 0.05; + dv[i] = 0.1; + } + if (powell(&rv, cal->devchan-1, dv, sr, 0.000001, 1000, ulimitfunc, + (void *)&cx, NULL, NULL) != 0) { + warning("icxUnderlyingLimit() failed for chan %d, ilimit %f\n",cal->devchan,ilimit); + return ilimit; + } + ulimitfunc((void *)&cx, dv); + + return cx.uilimit; +} + +/* ------------------------------------------------------ */ +/* Conversion and deltaE formular that include partial */ +/* derivatives, for use within fit parameter optimisations. */ + +/* CIE XYZ to perceptual Lab with partial derivatives. */ +void icxdXYZ2Lab(icmXYZNumber *w, double *out, double dout[3][3], double *in) { + double wp[3], tin[3], dtin[3]; + int i; + + wp[0] = w->X, wp[1] = w->Y, wp[2] = w->Z; + + for (i = 0; i < 3; i++) { + tin[i] = in[i]/wp[i]; + dtin[i] = 1.0/wp[i]; + + if (tin[i] > 0.008856451586) { + dtin[i] *= pow(tin[i], -2.0/3.0) / 3.0; + tin[i] = pow(tin[i],1.0/3.0); + } else { + dtin[i] *= 7.787036979; + tin[i] = 7.787036979 * tin[i] + 16.0/116.0; + } + } + + out[0] = 116.0 * tin[1] - 16.0; + dout[0][0] = 0.0; + dout[0][1] = 116.0 * dtin[1]; + dout[0][2] = 0.0; + + out[1] = 500.0 * (tin[0] - tin[1]); + dout[1][0] = 500.0 * dtin[0]; + dout[1][1] = 500.0 * -dtin[1]; + dout[1][2] = 0.0; + + out[2] = 200.0 * (tin[1] - tin[2]); + dout[2][0] = 0.0; + dout[2][1] = 200.0 * dtin[1]; + dout[2][2] = 200.0 * -dtin[2]; +} + + +/* Return the normal Delta E squared, given two Lab values, */ +/* including partial derivatives. */ +double icxdLabDEsq(double dout[2][3], double *Lab0, double *Lab1) { + double rv = 0.0, tt; + + tt = Lab0[0] - Lab1[0]; + dout[0][0] = 2.0 * tt; + dout[1][0] = -2.0 * tt; + rv += tt * tt; + tt = Lab0[1] - Lab1[1]; + dout[0][1] = 2.0 * tt; + dout[1][1] = -2.0 * tt; + rv += tt * tt; + tt = Lab0[2] - Lab1[2]; + dout[0][2] = 2.0 * tt; + dout[1][2] = -2.0 * tt; + rv += tt * tt; + return rv; +} + +/* Return the CIE94 Delta E color difference measure, squared */ +/* including partial derivatives. */ +double icxdCIE94sq(double dout[2][3], double Lab0[3], double Lab1[3]) { + double desq, _desq[2][3]; + double dlsq; + double dcsq, _dcsq[2][2]; /* == [x][1,2] */ + double c12, _c12[2][2]; /* == [x][1,2] */ + double dhsq, _dhsq[2][2]; /* == [x][1,2] */ + double rv; + + { + double dl, da, db; + + dl = Lab0[0] - Lab1[0]; + dlsq = dl * dl; /* dl squared */ + da = Lab0[1] - Lab1[1]; + db = Lab0[2] - Lab1[2]; + + + /* Compute normal Lab delta E squared */ + desq = dlsq + da * da + db * db; + _desq[0][0] = 2.0 * dl; + _desq[1][0] = -2.0 * dl; + _desq[0][1] = 2.0 * da; + _desq[1][1] = -2.0 * da; + _desq[0][2] = 2.0 * db; + _desq[1][2] = -2.0 * db; + } + + { + double c1, c2, dc, tt; + + /* Compute chromanance for the two colors */ + c1 = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]); + c2 = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]); + c12 = sqrt(c1 * c2); /* Symetric chromanance */ + + tt = 0.5 * (pow(c2, 0.5) + 1e-12)/(pow(c1, 1.5) + 1e-12); + _c12[0][0] = Lab0[1] * tt; + _c12[0][1] = Lab0[2] * tt; + tt = 0.5 * (pow(c1, 0.5) + 1e-12)/(pow(c2, 1.5) + 1e-12); + _c12[1][0] = Lab1[1] * tt; + _c12[1][1] = Lab1[2] * tt; + + /* delta chromanance squared */ + dc = c2 - c1; + dcsq = dc * dc; + if (c1 < 1e-12 || c2 < 1e-12) { + c1 += 1e-12; + c2 += 1e-12; + } + _dcsq[0][0] = -2.0 * Lab0[1] * (c2 - c1)/c1; + _dcsq[0][1] = -2.0 * Lab0[2] * (c2 - c1)/c1; + _dcsq[1][0] = 2.0 * Lab1[1] * (c2 - c1)/c2; + _dcsq[1][1] = 2.0 * Lab1[2] * (c2 - c1)/c2; + } + + /* Compute delta hue squared */ + dhsq = desq - dlsq - dcsq; + if (dhsq >= 0.0) { + _dhsq[0][0] = _desq[0][1] - _dcsq[0][0]; + _dhsq[0][1] = _desq[0][2] - _dcsq[0][1]; + _dhsq[1][0] = _desq[1][1] - _dcsq[1][0]; + _dhsq[1][1] = _desq[1][2] - _dcsq[1][1]; + } else { + dhsq = 0.0; + _dhsq[0][0] = 0.0; + _dhsq[0][1] = 0.0; + _dhsq[1][0] = 0.0; + _dhsq[1][1] = 0.0; + } + + { + double sc, scsq, scf; + double sh, shsq, shf; + + /* Weighting factors for delta chromanance & delta hue */ + sc = 1.0 + 0.048 * c12; + scsq = sc * sc; + + sh = 1.0 + 0.014 * c12; + shsq = sh * sh; + + rv = dlsq + dcsq/scsq + dhsq/shsq; + + scf = 0.048 * -2.0 * dcsq/(scsq * sc); + shf = 0.014 * -2.0 * dhsq/(shsq * sh); + dout[0][0] = _desq[0][0]; + dout[0][1] = _dcsq[0][0]/scsq + _c12[0][0] * scf + + _dhsq[0][0]/shsq + _c12[0][0] * shf; + dout[0][2] = _dcsq[0][1]/scsq + _c12[0][1] * scf + + _dhsq[0][1]/shsq + _c12[0][1] * shf; + dout[1][0] = _desq[1][0]; + dout[1][1] = _dcsq[1][0]/scsq + _c12[1][0] * scf + + _dhsq[1][0]/shsq + _c12[1][0] * shf; + dout[1][2] = _dcsq[1][1]/scsq + _c12[1][1] * scf + + _dhsq[1][1]/shsq + _c12[1][1] * shf; + return rv; + } +} + +// ~~99 not sure if these are correct: + +/* Return the normal Delta E given two Lab values, */ +/* including partial derivatives. */ +double icxdLabDE(double dout[2][3], double *Lab0, double *Lab1) { + double rv = 0.0, tt; + + tt = Lab0[0] - Lab1[0]; + dout[0][0] = 1.0 * tt; + dout[1][0] = -1.0 * tt; + rv += tt * tt; + tt = Lab0[1] - Lab1[1]; + dout[0][1] = 1.0 * tt; + dout[1][1] = -1.0 * tt; + rv += tt * tt; + tt = Lab0[2] - Lab1[2]; + dout[0][2] = 1.0 * tt; + dout[1][2] = -1.0 * tt; + rv += tt * tt; + return sqrt(rv); +} + +/* Return the CIE94 Delta E color difference measure */ +/* including partial derivatives. */ +double icxdCIE94(double dout[2][3], double Lab0[3], double Lab1[3]) { + double desq, _desq[2][3]; + double dlsq; + double dcsq, _dcsq[2][2]; /* == [x][1,2] */ + double c12, _c12[2][2]; /* == [x][1,2] */ + double dhsq, _dhsq[2][2]; /* == [x][1,2] */ + double rv; + + { + double dl, da, db; + + dl = Lab0[0] - Lab1[0]; + dlsq = dl * dl; /* dl squared */ + da = Lab0[1] - Lab1[1]; + db = Lab0[2] - Lab1[2]; + + + /* Compute normal Lab delta E squared */ + desq = dlsq + da * da + db * db; + _desq[0][0] = 1.0 * dl; + _desq[1][0] = -1.0 * dl; + _desq[0][1] = 1.0 * da; + _desq[1][1] = -1.0 * da; + _desq[0][2] = 1.0 * db; + _desq[1][2] = -1.0 * db; + } + + { + double c1, c2, dc, tt; + + /* Compute chromanance for the two colors */ + c1 = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]); + c2 = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]); + c12 = sqrt(c1 * c2); /* Symetric chromanance */ + + tt = 0.5 * (pow(c2, 0.5) + 1e-12)/(pow(c1, 1.5) + 1e-12); + _c12[0][0] = Lab0[1] * tt; + _c12[0][1] = Lab0[2] * tt; + tt = 0.5 * (pow(c1, 0.5) + 1e-12)/(pow(c2, 1.5) + 1e-12); + _c12[1][0] = Lab1[1] * tt; + _c12[1][1] = Lab1[2] * tt; + + /* delta chromanance squared */ + dc = c2 - c1; + dcsq = dc * dc; + if (c1 < 1e-12 || c2 < 1e-12) { + c1 += 1e-12; + c2 += 1e-12; + } + _dcsq[0][0] = -1.0 * Lab0[1] * (c2 - c1)/c1; + _dcsq[0][1] = -1.0 * Lab0[2] * (c2 - c1)/c1; + _dcsq[1][0] = 1.0 * Lab1[1] * (c2 - c1)/c2; + _dcsq[1][1] = 1.0 * Lab1[2] * (c2 - c1)/c2; + } + + /* Compute delta hue squared */ + dhsq = desq - dlsq - dcsq; + if (dhsq >= 0.0) { + _dhsq[0][0] = _desq[0][1] - _dcsq[0][0]; + _dhsq[0][1] = _desq[0][2] - _dcsq[0][1]; + _dhsq[1][0] = _desq[1][1] - _dcsq[1][0]; + _dhsq[1][1] = _desq[1][2] - _dcsq[1][1]; + } else { + dhsq = 0.0; + _dhsq[0][0] = 0.0; + _dhsq[0][1] = 0.0; + _dhsq[1][0] = 0.0; + _dhsq[1][1] = 0.0; + } + + { + double sc, scsq, scf; + double sh, shsq, shf; + + /* Weighting factors for delta chromanance & delta hue */ + sc = 1.0 + 0.048 * c12; + scsq = sc * sc; + + sh = 1.0 + 0.014 * c12; + shsq = sh * sh; + + rv = dlsq + dcsq/scsq + dhsq/shsq; + + scf = 0.048 * -1.0 * dcsq/(scsq * sc); + shf = 0.014 * -1.0 * dhsq/(shsq * sh); + dout[0][0] = _desq[0][0]; + dout[0][1] = _dcsq[0][0]/scsq + _c12[0][0] * scf + + _dhsq[0][0]/shsq + _c12[0][0] * shf; + dout[0][2] = _dcsq[0][1]/scsq + _c12[0][1] * scf + + _dhsq[0][1]/shsq + _c12[0][1] * shf; + dout[1][0] = _desq[1][0]; + dout[1][1] = _dcsq[1][0]/scsq + _c12[1][0] * scf + + _dhsq[1][0]/shsq + _c12[1][0] * shf; + dout[1][2] = _dcsq[1][1]/scsq + _c12[1][1] * scf + + _dhsq[1][1]/shsq + _c12[1][1] * shf; + return sqrt(rv); + } +} + +/* ------------------------------------------------------ */ +/* A power-like function, based on Graphics Gems adjustment curve. */ +/* Avoids "toe" problem of pure power. */ +/* Adjusted so that "power" 2 and 0.5 agree with real power at 0.5 */ + +double icx_powlike(double vv, double pp) { + double tt, g; + + if (pp >= 1.0) { + g = 2.0 * (pp - 1.0); + vv = vv/(g - g * vv + 1.0); + } else { + g = 2.0 - 2.0/pp; + vv = (vv - g * vv)/(1.0 - g * vv); + } + + return vv; +} + +/* Compute the necessary aproximate power, to transform */ +/* the given value from src to dst. They are assumed to be */ +/* in the range 0.0 .. 1.0 */ +double icx_powlike_needed(double src, double dst) { + double pp, g; + + if (dst <= src) { + g = -((src - dst)/(dst * src - dst)); + pp = (0.5 * g) + 1.0; + } else { + g = -((src - dst)/((dst - 1.0) * src)); + pp = 1.0/(1.0 - 0.5 * g); + } + + return pp; +} + +/* ------------------------------------------------------ */ +/* Parameterized transfer/dot gain function. */ +/* Used for device modelling. Including partial */ +/* derivative for input and parameters. */ + +/* NOTE that clamping the input values seems to cause */ +/* conjgrad() problems. */ + +/* Transfer function */ +double icxTransFunc( +double *v, /* Pointer to first parameter */ +int luord, /* Number of parameters */ +double vv /* Source of value */ +) { + double g; + int ord; + + /* Process all the shaper orders from low to high. */ + /* [These shapers were inspired by a Graphics Gem idea */ + /* (Gems IV, VI.3, "Fast Alternatives to Perlin's Bias and */ + /* Gain Functions, pp 401). */ + /* They have the nice properties that they are smooth, and */ + /* are monotonic. The control parameter has been */ + /* altered to have a range from -oo to +oo rather than 0.0 to 1.0 */ + /* so that the search space is less non-linear. */ + for (ord = 0; ord < luord; ord++) { + int nsec; /* Number of sections */ + double sec; /* Section */ + + g = v[ord]; /* Parameter */ + + nsec = ord + 1; /* Increase sections for each order */ + + vv *= (double)nsec; + + sec = floor(vv); + if (((int)sec) & 1) + g = -g; /* Alternate action in each section */ + vv -= sec; + if (g >= 0.0) { + vv = vv/(g - g * vv + 1.0); + } else { + vv = (vv - g * vv)/(1.0 - g * vv); + } + vv += sec; + vv /= (double)nsec; + } + + return vv; +} + +/* Inverse transfer function */ +double icxInvTransFunc( +double *v, /* Pointer to first parameter */ +int luord, /* Number of parameters */ +double vv /* Source of value */ +) { + double g; + int ord; + + /* Process the shaper orders in reverse from high to low. */ + /* [These shapers were inspired by a Graphics Gem idea */ + /* (Gems IV, VI.3, "Fast Alternatives to Perlin's Bias and */ + /* Gain Functions, pp 401). */ + /* They have the nice properties that they are smooth, and */ + /* are monotonic. The control parameter has been */ + /* altered to have a range from -oo to +oo rather than 0.0 to 1.0 */ + /* so that the search space is less non-linear. */ + for (ord = luord-1; ord >= 0; ord--) { + int nsec; /* Number of sections */ + double sec; /* Section */ + + g = -v[ord]; /* Inverse parameter */ + + nsec = ord + 1; /* Increase sections for each order */ + + vv *= (double)nsec; + + sec = floor(vv); + if (((int)sec) & 1) + g = -g; /* Alternate action in each section */ + vv -= sec; + if (g >= 0.0) { + vv = vv/(g - g * vv + 1.0); + } else { + vv = (vv - g * vv)/(1.0 - g * vv); + } + vv += sec; + vv /= (double)nsec; + } + + return vv; +} + +/* Transfer function with offset and scale */ +double icxSTransFunc( +double *v, /* Pointer to first parameter */ +int luord, /* Number of parameters */ +double vv, /* Source of value */ +double min, /* Scale values */ +double max +) { + max -= min; + + vv = (vv - min)/max; + vv = icxTransFunc(v, luord, vv); + vv = (vv * max) + min; + return vv; +} + +/* Inverse Transfer function with offset and scale */ +double icxInvSTransFunc( +double *v, /* Pointer to first parameter */ +int luord, /* Number of parameters */ +double vv, /* Source of value */ +double min, /* Scale values */ +double max +) { + max -= min; + + vv = (vv - min)/max; + vv = icxInvTransFunc(v, luord, vv); + vv = (vv * max) + min; + return vv; +} + +/* Transfer function with partial derivative */ +/* with respect to the parameters. */ +double icxdpTransFunc( +double *v, /* Pointer to first parameter */ +double *dv, /* Return derivative wrt each parameter */ +int luord, /* Number of parameters */ +double vv /* Source of value */ +) { + double g; + int i, ord; + + /* Process all the shaper orders from high to low. */ + for (ord = 0; ord < luord; ord++) { + double dsv; /* del for del in g */ + double ddv; /* del for del in vv */ + int nsec; /* Number of sections */ + double sec; /* Section */ + + g = v[ord]; /* Parameter */ + + nsec = ord + 1; /* Increase sections for each order */ + + vv *= (double)nsec; + + sec = floor(vv); + if (((int)sec) & 1) { + g = -g; /* Alternate action in each section */ + } + vv -= sec; + if (g >= 0.0) { + double tt = g - g * vv + 1.0; + dsv = (vv * vv - vv)/(tt * tt); + ddv = (g + 1.0)/(tt * tt); + vv = vv/tt; + } else { + double tt = 1.0 - g * vv; + dsv = (vv * vv - vv)/(tt * tt); + ddv = (1.0 - g)/(tt * tt); + vv = (vv - g * vv)/tt; + } + + vv += sec; + vv /= (double)nsec; + dsv /= (double)nsec; + if (((int)sec) & 1) + dsv = -dsv; + + dv[ord] = dsv; + for (i = ord - 1; i >= 0; i--) + dv[i] *= ddv; + } + + return vv; +} + +/* Transfer function with offset and scale, and */ +/* partial derivative with respect to the parameters. */ +double icxdpSTransFunc( +double *v, /* Pointer to first parameter */ +double *dv, /* Return derivative wrt each parameter */ +int luord, /* Number of parameters */ +double vv, /* Source of value */ +double min, /* Scale values */ +double max +) { + int i; + max -= min; + + vv = (vv - min)/max; + vv = icxdpTransFunc(v, dv, luord, vv); + vv = (vv * max) + min; + for (i = 0; i < luord; i++) + dv[i] *= max; + return vv; +} + + +/* Transfer function with partial derivative */ +/* with respect to the input value. */ +double icxdiTransFunc( +double *v, /* Pointer to first parameter */ +double *pdin, /* Return derivative wrt source value */ +int luord, /* Number of parameters */ +double vv /* Source of value */ +) { + double g, din; + int ord; + +#ifdef NEVER + if (vv < 0.0 || vv > 1.0) { + if (vv < 0.0) + vv = 0.0; + else + vv = 1.0; + + *pdin = 0.0; + return vv; + } +#endif + din = 1.0; + + /* Process all the shaper orders from high to low. */ + for (ord = 0; ord < luord; ord++) { + double ddv; /* del for del in vv */ + int nsec; /* Number of sections */ + double sec; /* Section */ + + g = v[ord]; /* Parameter */ + + nsec = ord + 1; /* Increase sections for each order */ + + vv *= (double)nsec; + + sec = floor(vv); + if (((int)sec) & 1) { + g = -g; /* Alternate action in each section */ + } + vv -= sec; + if (g >= 0.0) { + double tt = g - g * vv + 1.0; + ddv = (g + 1.0)/(tt * tt); + vv = vv/tt; + } else { + double tt = 1.0 - g * vv; + ddv = (1.0 - g)/(tt * tt); + vv = (vv - g * vv)/tt; + } + + vv += sec; + vv /= (double)nsec; + din *= ddv; + } + + *pdin = din; + return vv; +} + +/* Transfer function with offset and scale, and */ +/* partial derivative with respect to the input value. */ +double icxdiSTransFunc( +double *v, /* Pointer to first parameter */ +double *pdv, /* Return derivative wrt source value */ +int luord, /* Number of parameters */ +double vv, /* Source of value */ +double min, /* Scale values */ +double max +) { + max -= min; + + vv = (vv - min)/max; + vv = icxdiTransFunc(v, pdv, luord, vv); + vv = (vv * max) + min; + return vv; +} + + +/* Transfer function with partial derivative */ +/* with respect to the parameters and the input value. */ +double icxdpdiTransFunc( +double *v, /* Pointer to first parameter */ +double *dv, /* Return derivative wrt each parameter */ +double *pdin, /* Return derivative wrt source value */ +int luord, /* Number of parameters */ +double vv /* Source of value */ +) { + double g, din; + int i, ord; + +#ifdef NEVER + if (vv < 0.0 || vv > 1.0) { + if (vv < 0.0) + vv = 0.0; + else + vv = 1.0; + + for (ord = 0; ord < luord; ord++) + dv[ord] = 0.0; + *pdin = 0.0; + return vv; + } +#endif + din = 1.0; + + /* Process all the shaper orders from high to low. */ + for (ord = 0; ord < luord; ord++) { + double dsv; /* del for del in g */ + double ddv; /* del for del in vv */ + int nsec; /* Number of sections */ + double sec; /* Section */ + + g = v[ord]; /* Parameter */ + + nsec = ord + 1; /* Increase sections for each order */ + + vv *= (double)nsec; + + sec = floor(vv); + if (((int)sec) & 1) { + g = -g; /* Alternate action in each section */ + } + vv -= sec; + if (g >= 0.0) { + double tt = g - g * vv + 1.0; + dsv = (vv * vv - vv)/(tt * tt); + ddv = (g + 1.0)/(tt * tt); + vv = vv/tt; + } else { + double tt = 1.0 - g * vv; + dsv = (vv * vv - vv)/(tt * tt); + ddv = (1.0 - g)/(tt * tt); + vv = (vv - g * vv)/tt; + } + + vv += sec; + vv /= (double)nsec; + dsv /= (double)nsec; + if (((int)sec) & 1) + dsv = -dsv; + + dv[ord] = dsv; + for (i = ord - 1; i >= 0; i--) + dv[i] *= ddv; + din *= ddv; + } + + *pdin = din; + return vv; +} + +/* Transfer function with offset and scale, and */ +/* partial derivative with respect to the */ +/* parameters and the input value. */ +double icxdpdiSTransFunc( +double *v, /* Pointer to first parameter */ +double *dv, /* Return derivative wrt each parameter */ +double *pdin, /* Return derivative wrt source value */ +int luord, /* Number of parameters */ +double vv, /* Source of value */ +double min, /* Scale values */ +double max +) { + int i; + max -= min; + + vv = (vv - min)/max; + vv = icxdpdiTransFunc(v, dv, pdin, luord, vv); + vv = (vv * max) + min; + for (i = 0; i < luord; i++) + dv[i] *= max; + return vv; +} + +/* ------------------------------------------------------ */ +/* Multi-plane interpolation, used for device modelling. */ +/* Including partial derivative for input and parameters. */ +/* A simple flat plane is used for each output. */ + +/* Multi-plane interpolation - uses base + di slope values. */ +/* Parameters are assumed to be fdi groups of di + 1 parameters. */ +void icxPlaneInterp( +double *v, /* Pointer to first parameter [fdi * (di + 1)] */ +int fdi, /* Number of output channels */ +int di, /* Number of input channels */ +double *out, /* Resulting fdi values */ +double *in /* Input di values */ +) { + int e, f; + + for (f = 0; f < fdi; f++) { + for (out[f] = 0.0, e = 0; e < di; e++, v++) { + out[f] += in[e] * *v; + } + out[f] += *v; + } +} + + +/* Multii-plane interpolation with partial derivative */ +/* with respect to the input and parameters. */ +void icxdpdiPlaneInterp( +double *v, /* Pointer to first parameter value [fdi * (di + 1)] */ +double *dv, /* Return [1 + di] deriv. wrt each parameter v */ +double *din, /* Return [fdi * di] deriv. wrt each input value */ +int fdi, /* Number of output channels */ +int di, /* Number of input channels */ +double *out, /* Resulting fdi values */ +double *in /* Input di values */ +) { + int e, ee, f, g; + int dip2 = (di + 1); /* Output dim increment through parameters */ + + /* Compute the output values */ + for (f = 0; f < fdi; f++) { + for (out[f] = 0.0, e = 0; e < di; e++) + out[f] += in[e] * v[f * dip2 + e]; + out[f] += v[f * dip2 + e]; + } + + /* Since interpolation is verys simple, derivative are also simple */ + + /* Copy del for parameter to return array */ + for (e = 0; e < di; e++) + dv[e] = in[e]; + dv[e] = 1.0; + + /* Compute del of out[] from in[] */ + for (f = 0; f < fdi; f++) { + for (e = 0; e < di; e++) { + din[f * di + e] = v[f * dip2 + e]; + } + } +} + +/* ------------------------------------------------------ */ +/* Matrix cube interpolation, used for device modelling. */ +/* Including partial derivative for input and parameters. */ + +/* Matrix cube interpolation - interpolate between 2^di output corner values. */ +/* Parameters are assumed to be fdi groups of 2^di parameters. */ +void icxCubeInterp( +double *v, /* Pointer to first parameter */ +int fdi, /* Number of output channels */ +int di, /* Number of input channels */ +double *out, /* Resulting fdi values */ +double *in /* Input di values */ +) { + int e, f, g; + double gw[1 << MXDI]; /* weight for each matrix grid cube corner */ + + /* Compute corner weights needed for interpolation */ + gw[0] = 1.0; + for (e = 0, g = 1; e < di; e++, g *= 2) { + int i; + for (i = 0; i < g; i++) { + gw[g+i] = gw[i] * in[e]; + gw[i] *= (1.0 - in[e]); + } + } + + /* Now compute the output values */ + for (f = 0; f < fdi; f++) { + out[f] = 0.0; /* For each output value */ + for (e = 0; e < (1 << di); e++) { /* For all corners of cube */ + out[f] += gw[e] * *v; + v++; + } + } +} + + +/* Matrix cube interpolation. with partial derivative */ +/* with respect to the input and parameters. */ +void icxdpdiCubeInterp( +double *v, /* Pointer to first parameter value [fdi * 2^di] */ +double *dv, /* Return [2^di] deriv. wrt each parameter v */ +double *din, /* Return [fdi * di] deriv. wrt each input value */ +int fdi, /* Number of output channels */ +int di, /* Number of input channels */ +double *out, /* Resulting fdi values */ +double *in /* Input di values */ +) { + int e, ee, f, g; + int dip2 = (1 << di); + double gw[1 << MXDI]; /* weight for each matrix grid cube corner */ + + /* Compute corner weights needed for interpolation */ + gw[0] = 1.0; + for (e = 0, g = 1; e < di; e++, g *= 2) { + int i; + for (i = 0; i < g; i++) { + gw[g+i] = gw[i] * in[e]; + gw[i] *= (1.0 - in[e]); + } + } + + /* Now compute the output values */ + for (f = 0; f < fdi; f++) { + out[f] = 0.0; /* For each output value */ + for (ee = 0; ee < dip2; ee++) { /* For all corners of cube */ + out[f] += gw[ee] * v[f * dip2 + ee]; + } + } + + /* Copy del for parameter to return array */ + for (ee = 0; ee < dip2; ee++) { /* For all other corners of cube */ + dv[ee] = gw[ee]; /* del from parameter */ + } + + /* Compute del from in[] value we want */ + for (e = 0; e < di; e++) { /* For input we want del wrt */ + + for (f = 0; f < fdi; f++) + din[f * di + e] = 0.0; /* Zero del ready of accumulation */ + + for (ee = 0; ee < dip2; ee++) { /* For all corners of cube weights, */ + int e2; /* accumulate del from in[] we want. */ + double vv = 1.0; + + /* Compute in[] weighted cube corners for all except del of in[] we want */ + for (e2 = 0; e2 < di; e2++) { /* const from non del inputs */ + if (e2 == e) + continue; + if (ee & (1 << e2)) + vv *= in[e2]; + else + vv *= (1.0 - in[e2]); + } + + /* Accumulate contribution of in[] we want for corner to out[] we want */ + if (ee & (1 << e)) { + for (f = 0; f < fdi; f++) + din[f * di + e] += v[f * dip2 + ee] * vv; + } else { + for (f = 0; f < fdi; f++) + din[f * di + e] -= v[f * dip2 + ee] * vv; + } + } + } +} + +/* ------------------------------------------------------ */ +/* Matrix cube simplex interpolation, used for device modelling. */ + +/* Matrix cube simplex interpolation - interpolate between 2^di output corner values. */ +/* Parameters are assumed to be fdi groups of 2^di parameters. */ +void icxCubeSxInterp( +double *v, /* Pointer to first parameter */ +int fdi, /* Number of output channels */ +int di, /* Number of input channels */ +double *out, /* Resulting fdi values */ +double *in /* Input di values */ +) { + int si[MAX_CHAN]; /* in[] Sort index, [0] = smalest */ + +//{ +// double tout[MXDO]; +// +// icxCubeInterp(v, fdi, di, tout, in); +//printf("\n~1 Cube interp result = %f\n",tout[0]); +//} + +//printf("~1 icxCubeSxInterp: %f %f %f\n", in[0], in[1], in[2]); + /* Do insertion sort on coordinates, smallest to largest. */ + { + int ff, vf; + unsigned int e; + double v; + for (e = 0; e < di; e++) + si[e] = e; /* Initial unsorted indexes */ + + for (e = 1; e < di; e++) { + ff = e; + v = in[si[ff]]; + vf = ff; + while (ff > 0 && in[si[ff-1]] > v) { + si[ff] = si[ff-1]; + ff--; + } + si[ff] = vf; + } + } +//printf("~1 sort order %d %d %d\n", si[0], si[1], si[2]); +//printf(" from %f %f %f\n", in[si[0]], in[si[1]], in[si[2]]); + + /* Now compute the weightings, simplex vertices and output values */ + { + unsigned int e, f; + double w; /* Current vertex weight */ + + w = 1.0 - in[si[di-1]]; /* Vertex at base of cell */ + for (f = 0; f < fdi; f++) { + out[f] = w * v[f * (1 << di)]; +//printf("~1 out[%d] = %f = %f * %f\n",f,out[f],w,v[f * (1 << di)]); + } + + for (e = di-1; e > 0; e--) { /* Middle verticies */ + w = in[si[e]] - in[si[e-1]]; + v += (1 << si[e]); /* Move to top of cell in next largest dimension */ + for (f = 0; f < fdi; f++) { + out[f] += w * v[f * (1 << di)]; +//printf("~1 out[%d] = %f += %f * %f\n",f,out[f],w,v[f * (1 << di)]); + } + } + + w = in[si[0]]; + v += (1 << si[0]); /* Far corner from base of cell */ + for (f = 0; f < fdi; f++) { + out[f] += w * v[f * (1 << di)]; +//printf("~1 out[%d] = %f += %f * %f\n",f,out[f],w,v[f * (1 << di)]); + } + } +} + +/* ------------------------------------------------------ */ +/* Matrix multiplication, used for device modelling. */ +/* Including partial derivative for input and parameters. */ + + +/* 3x3 matrix multiplication, with the matrix in a 1D array */ +/* with respect to the input and parameters. */ +void icxMulBy3x3Parm( + double out[3], /* Return input multiplied by matrix */ + double mat[9], /* Matrix organised in [slow][fast] order */ + double in[3] /* Input values */ +) { + double *v = mat, ov[3]; + int e, f; + + /* Compute the output values */ + for (f = 0; f < 3; f++) { + ov[f] = 0.0; /* For each output value */ + for (e = 0; e < 3; e++) { + ov[f] += *v++ * in[e]; + } + } + out[0] = ov[0]; + out[1] = ov[1]; + out[2] = ov[2]; +} + + +/* 3x3 matrix multiplication, with partial derivatives */ +/* with respect to the input and parameters. */ +void icxdpdiMulBy3x3Parm( + double out[3], /* Return input multiplied by matrix */ + double dv[3][9], /* Return deriv for each [output] with respect to [param] */ + double din[3][3], /* Return deriv for each [output] with respect to [input] */ + double mat[9], /* Matrix organised in [slow][fast] order */ + double in[3] /* Input values */ +) { + double *v, ov[3]; + int e, f; + + /* Compute the output values */ + v = mat; + for (f = 0; f < 3; f++) { + ov[f] = 0.0; /* For each output value */ + for (e = 0; e < 3; e++) { + ov[f] += *v++ * in[e]; + } + } + + /* Compute deriv. with respect to the matrix parameter % 3 */ + /* This is pretty simple for a matrix ... */ + for (f = 0; f < 3; f++) { + for (e = 0; e < 9; e++) { + if (e/3 == f) + dv[f][e] = in[e % 3]; + else + dv[f][e] = 0.0; + } + } + + /* Compute deriv. with respect to the input values */ + /* This is pretty simple for a matrix ... */ + v = mat; + for (f = 0; f < 3; f++) + for (e = 0; e < 3; e++) + din[f][e] = *v++; + + out[0] = ov[0]; + out[1] = ov[1]; + out[2] = ov[2]; +} + +#undef stricmp + + + + + + + + + + + + + + + + + + + + + + + |