diff options
Diffstat (limited to 'gamut/nearsmth.c')
| -rw-r--r-- | gamut/nearsmth.c | 3570 |
1 files changed, 3570 insertions, 0 deletions
diff --git a/gamut/nearsmth.c b/gamut/nearsmth.c new file mode 100644 index 0000000..1f48f20 --- /dev/null +++ b/gamut/nearsmth.c @@ -0,0 +1,3570 @@ + +/* + * nearsmth + * + * Gamut mapping support routine that creates a list of + * guide vectors that map from a source to destination + * gamut, smoothed to retain reasonably even spacing. + * + * Author: Graeme W. Gill + * Date: 17/1/2002 + * Version: 1.00 + * + * Copyright 2002 - 2006 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. + */ + +/* + Description: + + We create a set of "guide vectors" that map the source gamut to + the destination, for use by the gammap code in creating + a 3D gamut mapping. + + (See gammap.txt for a more detailed descrition) + + */ + +/* + * TTBD: + * + * It might work better if the cusp mapping had separate control + * over the L and h degree of map, as well as the L and h effective radius ? + * That way, saturation hue distortions with L might be reduced. + * + * Improve error handling. + * + * Major defect with some gamut combinations is "button" around + * cusps. Not sure what the mechanism is, since it's not obvious + * from the 3D vector plots what the cause is. (fixed ?) + * Due to poor internal control ? + * + * Mapping to very small, odd shaped gamuts (ie. Bonet) is poor - + * there are various bugs and artefacts to be figured out. + */ + +#include <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +#include <fcntl.h> +#include <string.h> +#include <math.h> +#include "icc.h" +#include "numlib.h" +#include "rspl.h" +#include "gamut.h" +#include "nearsmth.h" +#include "vrml.h" + +#undef SAVE_VRMLS /* [Und] Save various vrml's */ +#undef PLOT_MAPPING_INFLUENCE /* [Und] Plot sci_gam colored by dominant guide influence: */ + /* Absolute = red, Relative = yellow, Radial = blue, Depth = green */ +#undef PLOT_AXES /* [Und] */ +#undef PLOT_EVECTS /* [Und] Create VRML of error correction vectors */ +#undef VERB /* [Und] [0] If <= 1, print progress headings */ + /* if > 1, print information about everything */ +#undef SHOW_NEIGB_WEIGHTS /* [Und] Show the weighting for each point of neighbours */ + +#undef DIAG_POINTS /* [Und] Short circuite mapping and show vectors of various */ + /* intermediate points (see #ifdef DIAG_POINTS) */ + +#undef PLOT_DIGAM /* [Und] Rather than DST_GMT - don't free it (#def in gammap.c too) */ + +#define SUM_POW 2.0 /* Delta's are sum of component deltas ^ SUM_POW */ +#define LIGHT_L 70.0 /* "light" L/J value */ +#define DARK_L 5.0 /* "dark" L/J value */ +#define NEUTRAL_C 20.0 /* "neutral" C value */ +#define NO_TRIALS 6 /* [6] Number of random trials */ +#define VECSMOOTHING /* [Def] Enable vector smoothing */ +#define VECADJPASSES 3 /* [3] Adjust vectors after smoothing to be on dest gamut */ +#define RSPLPASSES 4 /* [4] Number of rspl adjustment passes */ +#define RSPLSCALE 1.8 /* [1.8] Offset within gamut for rspl smoothingto aim for */ +#define SHRINK 5.0 /* Shrunk destination evect surface factor */ +#define CYLIN_SUBVEC /* [Def] Make sub-vectors always cylindrical direction */ +#define SUBVEC_SMOOTHING /* [Def] Smooth the sub-vectors */ + + /* Experimental - not used: */ + + /* This has similar effects to lowering SUM_POW without the side effects */ + /* and improves hue detail for small destination gamuts. */ + /* (This and lxpow are pretty hacky. Is there a better way ?) */ +#undef EMPH_NEUTRAL //0.5 /* Emphasis strength near neutral */ +#undef EMPH_THR //10.0 /* delta C threshold above which it kicks in */ + +#undef LINEAR_HUE_SUM /* Make delta^2 = (sqrt(l^2 + c^2) + h)^2 */ + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +#if defined(VERB) +# define VA(xxxx) printf xxxx +# if VERB > 1 +# define VB(xxxx) printf xxxx +# else +# define VB(xxxx) +# endif +#else +# define VA(xxxx) +# define VB(xxxx) +#endif + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +#if defined(SAVE_VRMLS) && defined(PLOT_MAPPING_INFLUENCE) +static void create_influence_plot(nearsmth *smp, int nmpts); +#endif + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Compute the weighted delta E squared of in1 - in2 */ +/* (This is like the CIE DE94) */ +static double wdesq( +double in1[3], /* Destination location */ +double in2[3], /* Source location */ +double lweight, +double cweight, +double hweight, +double sumpow /* Sum power. 0.0 == 2.0 */ +) { + double desq, dhsq; + double dlsq, dcsq; + double vv; + double dc, c1, c2; + +//printf("~1 wdesq got %f %f %f and %f %f %f\n", in1[0], in1[1], in1[2], in2[0], in2[1], in2[2]); + /* Compute delta L squared and delta E squared */ + { + double dl, da, db; + dl = in1[0] - in2[0]; + dlsq = dl * dl; /* dl squared */ + da = in1[1] - in2[1]; + db = in1[2] - in2[2]; + + desq = dlsq + da * da + db * db; + } + + /* compute delta chromanance squared */ + { + + /* Compute chromanance for the two colors */ + c1 = sqrt(in1[1] * in1[1] + in1[2] * in1[2]); + c2 = sqrt(in2[1] * in2[1] + in2[2] * in2[2]); + + dc = c1 - c2; + dcsq = dc * dc; + + /* [ Making dcsq = sqrt(dcsq) here seemes */ + /* to improve the saturation result. Subsumed by a.xl ? ] */ + } + + /* Compute delta hue squared */ + if ((dhsq = desq - dlsq - dcsq) < 0.0) + dhsq = 0.0; + +#ifdef EMPH_NEUTRAL /* Emphasise hue differences whenc dc is large and we are */ + /* close to the neutral axis */ + vv = 3.0 / (1.0 + 0.03 * c1); /* Full strength scale factor from dest location */ + vv = 1.0 + EMPH_NEUTRAL * (vv - 1.0); /* Reduced strength scale factor */ + vv *= (dc + EMPH_THR)/EMPH_THR; + dhsq *= vv * vv; /* Scale squared hue delta */ +#endif + + if (sumpow == 0.0 || sumpow == 2.0) { /* Normal sum of squares */ +#ifdef HACK + vv = sqrt(lweight * dlsq + cweight * dcsq) + sqrt(hweight * dhsq); + vv *= vv; + vv = fabs(vv); /* Avoid -0.0 */ +#else + vv = lweight * dlsq + cweight * dcsq + hweight * dhsq; + vv = fabs(vv); /* Avoid -0.0 */ +#endif + } else { + sumpow *= 0.5; + vv = lweight * pow(dlsq, sumpow) + cweight * pow(dcsq,sumpow) + hweight * pow(dhsq,sumpow); + vv = fabs(vv); /* Avoid -0.0 */ + vv = pow(vv, 1.0/sumpow); + } + +//printf("~1 returning wdesq %f from %f * %f + %f * %f + %f * %f\n", fabs(vv),lweight, dlsq, cweight, dcsq, hweight, dhsq); + + return vv; +} + +/* Compute the LCh differences squared of in1 - in2 */ +/* (This is like the CIE DE94) */ +static void diffLCh( +double out[3], +double in1[3], /* Destination location */ +double in2[3] /* Source location */ +) { + double desq, dhsq; + double dlsq, dcsq; + double vv; + double dc, c1, c2; + + /* Compute delta L squared and delta E squared */ + { + double dl, da, db; + dl = in1[0] - in2[0]; + dlsq = dl * dl; /* dl squared */ + da = in1[1] - in2[1]; + db = in1[2] - in2[2]; + + desq = dlsq + da * da + db * db; + } + + /* compute delta chromanance squared */ + { + + /* Compute chromanance for the two colors */ + c1 = sqrt(in1[1] * in1[1] + in1[2] * in1[2]); + c2 = sqrt(in2[1] * in2[1] + in2[2] * in2[2]); + + dc = c1 - c2; + dcsq = dc * dc; + + /* [ Making dcsq = sqrt(dcsq) here seemes */ + /* to improve the saturation result. Subsumed by a.xl ? ] */ + } + + /* Compute delta hue squared */ + if ((dhsq = desq - dlsq - dcsq) < 0.0) + dhsq = 0.0; + +#ifdef EMPH_NEUTRAL /* Emphasise hue differences whenc dc is large and we are */ + /* close to the neutral axis */ + vv = 3.0 / (1.0 + 0.03 * c1); /* Full strength scale factor from dest location */ + vv = 1.0 + EMPH_NEUTRAL * (vv - 1.0); /* Reduced strength scale factor */ + vv *= (dc + EMPH_THR)/EMPH_THR; + dhsq *= vv * vv; /* Scale squared hue delta */ +#endif + + out[0] = dlsq; + out[1] = dcsq; + out[2] = dhsq; +} + +/* Given the weighting structure and the relevant point locations */ +/* return the total weighted error squared. */ +static double comperr( +nearsmth *p, /* Guide point data */ +gammapweights *w, /* weightings */ +double dtp[3], /* Dest test point being evaluated */ +double aodv[3], /* Weighted destination closest value to source */ +double drv[3], /* Source mapped radially to dest */ +double dcratio, /* Depth compression ratio of mapping */ +double dxratio /* Depth expansion ratio of mapping */ +) { + double a_o; + double va, vr = 0.0, vl, vd, vv = 0.0; + + /* Absolute, Delta E^2 between test point and destination closest */ + /* aodv is already positioned acording to the LCh weights, */ + /* so weight as per average of these */ + a_o = w->a.o; + va = wdesq(dtp, aodv, a_o, a_o, a_o, SUM_POW); + + /* Radial. Delta E^2 between test point and source mapped radially to dest gamut */ + vl = wdesq(dtp, drv, w->rl.l, w->rl.c, w->rl.h, SUM_POW); + + /* Depth ratio error^2. */ + vd = w->d.co * dcratio * dcratio + + w->d.xo * dxratio * dxratio; + + /* Diagnostic values */ + p->dbgv[0] = va; + p->dbgv[1] = vr; + p->dbgv[2] = vl; + p->dbgv[3] = vd; + + vv = va + vr + vl + vd; /* Sum of squares */ +// vv = sqrt(va) + sqrt(vr) + sqrt(vl) + sqrt(vd); /* Linear sum is better ? */ +// vv = pow(va, 0.7) + pow(vr, 0.7) + pow(vl, 0.7) + pow(vd, 0.7); /* Linear sum is better ? */ + +#ifdef NEVER + printf("~1 dtp %f %f %f\n", dtp[0], dtp[1], dtp[2]); + printf("~1 va = %f from aodv %f %f %f, weight %f\n", va, aodv[0], aodv[1], aodv[2], a_o); + printf("~1 vl = %f from drv %f %f %f, weights %f %f %f\n", vl, drv[0], drv[1], drv[2], w->rl.l, w->rl.c, w->rl.h); + printf("~1 vd = %f from d.co %f d.xo %f, weights %f %f\n", vd, w->d.co,w->d.xo,dcratio * dcratio,dxratio * dxratio); + printf("~1 return vv = %f\n", vv); +#endif /* NEVER */ + + return vv; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/* Structure to hold context for powell optimisation */ +/* and cusp mapping function. */ +struct _smthopt { + /* optimisation */ + int pass; /* Itteration round */ + int ix; /* Index of point being optimized */ + nearsmth *p; /* Point being optimised */ + int useexp; /* Flag indicating whether expansion is permitted */ + int debug; /* debug flag */ + gamut *shgam; /* for optfunc1a */ + + /* Setup state */ + int isJab; /* Flag indicating Jab rather than Lab space */ + int donaxis; /* Flag indicating whether neutral axis information is real */ + int docusp; /* Flag indicating whether cusp information is present */ + gammapweights *xwh; /* Structure holding expanded hextant weightings */ + gamut *sgam; /* Source colorspace gamut */ + + /* Cusp alignment mapping */ + /* 0 = src, 1 = dst, then cusp then value(s) */ + double cusps[2][9][3]; /* raw cusp values - red .. magenta, white [6], black [7] & grey [8] */ + double rot[2][3][4]; /* Rotation to align to black/white center */ + double irot[2][3][4]; /* Inverse rotation */ + double cusp_lab[2][9][3]; /* Cusp + B&W + grey rotated Lab value */ + double cusp_lch[2][6][3]; /* Cusp LCH value */ + double cusp_pe[2][6][4]; /* L direction plane equations per segment */ + double cusp_bc[2][6][2][3][3]; /* light/dark to/from 3x3 baricentic transform matrix */ + + /* Inversion support */ + double tv[3]; + gammapweights *wt; /* Weights for this inversion */ + + double mm[3][4]; /* Direction alignment rotation */ + double m2[2][2]; /* Additional matrix to alight a with L axis */ + double manv[3]; /* anv[] transformed by mm and m2 */ + +}; typedef struct _smthopt smthopt; + +static void init_ce(smthopt *s, gamut *sc_gam, gamut *d_gam, int src_kbp, int dst_kbp, double d_bp[3]); +static void comp_ce(smthopt *s, double out[3], double in[3], gammapweights *wt); +static void inv_comp_ce(smthopt *s, double out[3], double in[3], gammapweights *wt); +static double comp_naxbf(smthopt *s, double in[3]); +static double comp_lvc(smthopt *s, double in[3]); + +static double spow(double arg, double ex) { + if (arg < 0.0) + return -pow(-arg, ex); + else + return pow(arg, ex); +} + +static void spow3(double *out, double *in, double ex) { + int j; + for (j = 0; j < 3; j++) { + if (in[j] < 0.0) + out[j] = -pow(-in[j], ex); + else + out[j] = pow(in[j], ex); + } +} + +/* Powell optimisation function for setting minimal absolute error target point. */ +/* We get a 2D plane in the 3D space, of the destination point, */ +/* who's location we are optimizing. */ +static double optfunc1( +void *fdata, +double *_dv +) { + smthopt *s = (smthopt *)fdata; + nearsmth *p = s->p; /* Point being optimised */ + int i, j, k; + double dv[3]; /* 3D point in question */ + double ddv[3]; /* Point in question mapped to dst surface */ + double delch[3]; + double rv; /* Out of gamut, return value */ + + /* Convert from 2D to 3D. */ + dv[2] = _dv[1]; + dv[1] = _dv[0]; + dv[0] = 50.0; + icmMul3By3x4(dv, p->m3d, dv); + +//printf("~1 optfunc1 got 2D %f %f -> 3D %f %f %f\n", _dv[0], _dv[1], dv[0], dv[1], dv[2]); + p->dgam->radial(p->dgam, ddv, dv); /* Map to dst surface to check current location */ +//printf("~1 optfunc1 got %f %f %f -> surface %f %f %f\n", dv[0], dv[1], dv[2], ddv[0], ddv[1], ddv[2]); + + if (p->swap) { + /* This is actually a point on the real source gamut, so */ + /* convert to cusp mapped rotated, elevated source gamut value */ + comp_ce(s, ddv, ddv, &p->wt); +// printf("~1 after rot & elevate got %f %f %f\n",ddv[0],ddv[1],ddv[2]); + } + +#ifdef NEVER + /* Absolute weighted delta E between source and dest test point */ + rv = wdesq(ddv, p->sv, p->wt.ra.l, p->wt.ra.c, p->wt.ra.h, SUM_POW); +#else + { + double ppp = p->wt.a.lxpow; + double thr = p->wt.a.lxthr; /* Xover between normal and power */ + double sumpow = SUM_POW; + + diffLCh(delch, ddv, p->sv); + + if (sumpow == 0.0 || sumpow == 2.0) { /* Normal sum of squares */ +#ifdef LINEAR_HUE_SUM + double ll, cc, hh; + ll = p->wt.ra.l * pow(delch[0], ppp) * thr/pow(thr, ppp); + cc = p->wt.ra.c * delch[1]; + hh = p->wt.ra.h * delch[2]; + rv = sqrt(ll + cc) + sqrt(hh); + rv *= rv; +#else + rv = p->wt.ra.l * pow(delch[0], ppp) * thr/pow(thr, ppp) + + p->wt.ra.c * delch[1] + + p->wt.ra.h * delch[2]; +#endif + } else { + sumpow *= 0.5; + + rv = p->wt.ra.l * pow(delch[0], ppp * sumpow) * thr/pow(thr, ppp * sumpow) + + p->wt.ra.c * pow(delch[1], sumpow) + + p->wt.ra.h * pow(delch[2], sumpow); + } + } +#endif + + if (s->debug) + printf("debug: rv = %f from %f %f %f\n",rv, ddv[0], ddv[1], ddv[2]); + +//printf("~1 sv %4.2f %4.2f %4.2f, ddv %4.2f %4.2f %4.2f\n", p->sv[0], p->sv[1], p->sv[2], ddv[0], ddv[1], ddv[2]); +//printf("~1 rv = %f\n",rv); + return rv; +} + +/* Powell optimisation function for setting minimal absolute error target point, */ +/* from dest gamut to shrunk destination gamut. */ +/* We get a 2D plane in the 3D space, of the destination point, */ +/* who's location we are optimizing. */ +static double optfunc1a( +void *fdata, +double *_dv +) { + smthopt *s = (smthopt *)fdata; + nearsmth *p = s->p; /* Point being optimised */ + int i, j, k; + double dv[3]; /* 3D point in question */ + double ddv[3]; /* Point in question mapped to shgam surface */ + double delch[3]; + double rv; /* Out of gamut, return value */ + + /* Convert from 2D to 3D. */ + dv[2] = _dv[1]; + dv[1] = _dv[0]; + dv[0] = 50.0; + icmMul3By3x4(dv, p->m3d, dv); + +//printf("~1 optfunc1a got 2D %f %f -> 3D %f %f %f\n", _dv[0], _dv[1], dv[0], dv[1], dv[2]); + s->shgam->radial(s->shgam, ddv, dv); /* Map to shgam surface to check current location */ +//printf("~1 optfunc1a got %f %f %f -> surface %f %f %f\n", dv[0], dv[1], dv[2], ddv[0], ddv[1], ddv[2]); + +#ifdef NEVER + /* Absolute weighted delta E between source and dest test point */ + rv = wdesq(ddv, p->sv, p->wt.ra.l, p->wt.ra.c, p->wt.ra.h, SUM_POW); +#else + { + double ppp = p->wt.a.lxpow; + double thr = p->wt.a.lxthr; /* Xover between normal and power */ + double sumpow = SUM_POW; + + diffLCh(delch, ddv, p->dv); + + if (sumpow == 0.0 || sumpow == 2.0) { /* Normal sum of squares */ +#ifdef LINEAR_HUE_SUM + double ll, cc, hh; + ll = p->wt.ra.l * pow(delch[0], ppp) * thr/pow(thr, ppp); + cc = p->wt.ra.c * delch[1]; + hh = p->wt.ra.h * delch[2]; + rv = sqrt(ll + cc) + sqrt(hh); + rv *= rv; +#else + rv = p->wt.ra.l * pow(delch[0], ppp) * thr/pow(thr, ppp) + + p->wt.ra.c * delch[1] + + p->wt.ra.h * delch[2]; +#endif + } else { + sumpow *= 0.5; + + rv = p->wt.ra.l * pow(delch[0], ppp * sumpow) * thr/pow(thr, ppp * sumpow) + + p->wt.ra.c * pow(delch[1], sumpow) + + p->wt.ra.h * pow(delch[2], sumpow); + } + } +#endif + + if (s->debug) + printf("debug: rv = %f from %f %f %f\n",rv, ddv[0], ddv[1], ddv[2]); + +//printf("~1 sv %4.2f %4.2f %4.2f, ddv %4.2f %4.2f %4.2f\n", p->sv[0], p->sv[1], p->sv[2], ddv[0], ddv[1], ddv[2]); +//printf("~1 rv = %f\n",rv); + return rv; +} + + +/* Compute available depth errors p->dcratio and p->dxratio */ +static void comp_depth( +smthopt *s, +nearsmth *p, /* Point being optimized */ +double *dv /* 3D Location being evaluated */ +) { + double *sv, nv[3], nl; /* Source, dest points, normalized vector between them */ + double mint, maxt; + gtri *mintri = NULL, *maxtri = NULL; + + sv = p->_sv; + + p->dcratio = p->dxratio = 0.0; /* default, no depth error */ + + icmSub3(nv, dv, sv); /* Mapping vector */ + nl = icmNorm3(nv); /* It's length */ + if (nl > 0.1) { /* If mapping is non trivial */ + + icmScale3(nv, nv, 1.0/nl); /* Make mapping vector normal */ + + /* Compute actual depth of ray into destination (norm) or from source (expansion) gamut */ + if (p->dgam->vector_isect(p->dgam, sv, dv, NULL, NULL, &mint, &maxt, &mintri, &maxtri) != 0) { + double angle; + + /* The scale factor discounts the depth ratio as the mapping */ + /* vector gets more angled. It has a sin^2 characteristic */ + /* This is so that the depth error has some continuity if it */ + /* gets closer to being parallel to the destination gamut surface. */ + +//printf("\n~1 ix %d: %f %f %f -> %f %f %f\n isect at t %f and %f\n", s->ix, sv[0], sv[1], sv[2], dv[0], dv[1], dv[2], mint, maxt); + p->gflag = p->vflag = 0; + + if (mint < -1e-8 && maxt < -1e-8) { + p->gflag = 1; /* Gamut compression but */ + p->vflag = 2; /* vector is expanding */ + + } else if (mint > 1e-8 && maxt > -1e-8) { + p->gflag = 1; /* Gamut compression and */ + p->vflag = 1; /* vector compression */ + angle = icmDot3(nv, mintri->pe); + angle *= angle; /* sin squared */ + p->dcratio = angle * 2.0/(maxt + mint - 2.0); +//printf("~1 %d: comp depth ratio %f, angle %f\n", s->ix, p->dratio, angle); + + } else if (mint < -1e-8 && maxt > -1e-8) { + if (fabs(mint) < (fabs(maxt) - 1e-8)) { + p->gflag = 2; /* Gamut expansion but */ + p->vflag = 1; /* vector is compressing */ + + } else if (fabs(mint) > (fabs(maxt) + 1e-8)) { + p->gflag = 2; /* Gamut expansion and */ + p->vflag = 2; /* vector is expanding */ + angle = icmDot3(nv, maxtri->pe); + angle *= angle; /* sin squared */ + p->dxratio = angle * 2.0/-mint; +//printf("~1 %d: exp depth ratio %f, angle %f\n", s->ix, p->dratio, angle); + + } + } + } + } +} + +/* Powell optimisation function for non-relative error optimization. */ +/* We get a 2D point in the 3D space. */ +static double optfunc2( +void *fdata, +double *_dv +) { + smthopt *s = (smthopt *)fdata; + nearsmth *p = s->p; /* Point being optimised */ + double dv[3], ddv[3]; /* Dest point */ + double rv; /* Return value */ + + /* Convert from 2D to 3D. */ + dv[2] = _dv[1]; + dv[1] = _dv[0]; + dv[0] = 50.0; + icmMul3By3x4(dv, p->m3d, dv); +//printf("~1 optfunc2 got 2D %f %f -> 3D %f %f %f\n", _dv[0], _dv[1], dv[0], dv[1], dv[2]); + + p->dgam->radial(p->dgam, ddv, dv); /* Map to dst surface to check current location */ +//printf("~1 optfunc2 got %f %f %f -> surface %f %f %f\n", dv[0], dv[1], dv[2], ddv[0], ddv[1], ddv[2]); + +//printf("~1 optfunc2 sv %4.2f %4.2f %4.2f, dv %4.2f %4.2f %4.2f\n", p->sv[0], p->sv[1], p->sv[2], ddv[0], ddv[1], ddv[2]); + + /* Compute available depth errors p->dcratio and p->dxratio */ + comp_depth(s, p, ddv); + + /* Compute weighted delta E being minimised. */ + rv = comperr(p, &p->wt, ddv, p->aodv, p->drv, p->dcratio, p->dxratio); + + if (s->debug) { + printf("~1 sv = %f %f %f\n", p->sv[0], p->sv[1], p->sv[2]); + printf("~1 dv = %f %f %f\n", ddv[0], ddv[1], ddv[2]); + printf("~1 aodv = %f %f %f\n", p->aodv[0], p->aodv[1], p->aodv[2]); + printf("~1 drv = %f %f %f\n", p->drv[0], p->drv[1], p->drv[2]); + printf("debug:%d: rv = %f from %f %f %f\n",s->ix, rv, dv[0], dv[1], dv[2]); + } + +//printf("~1 rv = %f from %f %f\n",rv, _dv[0], _dv[1]); + +//printf("~1 rv = %f\n\n",rv); + return rv; +} + +/* -------------------------------------------- */ + +/* Setup the neutral axis and cusp mapping structure information */ +static void init_ce( +smthopt *s, /* Context for cusp mapping being set. */ +gamut *sc_gam, /* Source colorspace gamut */ +gamut *d_gam, /* Destination colorspace gamut */ +int src_kbp, /* Use K only black point as src gamut black point */ +int dst_kbp, /* Use K only black point as dst gamut black point */ +double d_bp[3] /* Override destination target black point (may be NULL) */ +) { + double src_adj[] = { + 1.1639766020018968e+224, 1.0605092189369252e-153, 3.5252483622572622e+257, + 1.3051549117649167e+214, 3.2984590678749676e-033, 1.8786244212510033e-153, + 1.2018790902224465e+049, 1.0618629743651763e-153, 5.5513445545255624e+233, + 3.3509081077514219e+242, 2.0076462988863408e-139, 3.2823498214286135e-318, + 7.7791723264448801e-260, 9.5956158769288055e+281, 2.5912667577703660e+161, + 5.2030128643503829e-085, 5.8235640814905865e+180, 4.0784546104861075e-033, + 3.6621812661291286e+098, 1.6417826055515754e-086, 8.2656018530749330e+097, + 9.3028116527073026e+242, 2.9127574654725916e+180, 1.9984697356129145e-139, + -2.1117351731638832e+003 }; + double saval; + int sd; + int i, j, k; + + VA(("init_ce called\n")); + + s->donaxis = 1; /* Assume real neutral axis info */ + s->docusp = 1; /* Assume real cusp info */ + + s->isJab = sc_gam->isJab; + + /* Set some default values for src white/black/grey */ + + /* Get the white and black point info */ + if (src_kbp) { + if (sc_gam->getwb(sc_gam, NULL, NULL, NULL, s->cusps[0][6], NULL, s->cusps[0][7]) != 0) { + VB(("getting src wb points failed\n")); + s->cusps[0][6][0] = 100.0; + s->cusps[0][7][0] = 0.0; + s->cusps[0][8][0] = 50.0; + s->donaxis = 0; + } + } else { + if (sc_gam->getwb(sc_gam, NULL, NULL, NULL, s->cusps[0][6], s->cusps[0][7], NULL) != 0) { + VB(("getting src wb points failed\n")); + s->cusps[0][6][0] = 100.0; + s->cusps[0][7][0] = 0.0; + s->cusps[0][8][0] = 50.0; + s->donaxis = 0; + } + } + + if (dst_kbp) { + if (d_gam->getwb(d_gam, NULL, NULL, NULL, s->cusps[1][6], NULL, s->cusps[1][7]) != 0) { + VB(("getting dest wb points failed\n")); + s->cusps[1][6][0] = 100.0; + s->cusps[1][7][0] = 0.0; + s->cusps[1][8][0] = 50.0; + s->donaxis = 0; + } + } else { + if (d_gam->getwb(d_gam, NULL, NULL, NULL, s->cusps[1][6], s->cusps[1][7], NULL) != 0) { + VB(("getting dest wb points failed\n")); + s->cusps[1][6][0] = 100.0; + s->cusps[1][7][0] = 0.0; + s->cusps[1][8][0] = 50.0; + s->donaxis = 0; + } + } + if (d_bp != NULL) { /* Use override destination black point */ + icmCpy3(s->cusps[1][7], d_bp); + } + + /* Get the cusp info */ + if (sc_gam->getcusps(sc_gam, s->cusps[0]) != 0 || d_gam->getcusps(d_gam, s->cusps[1]) != 0) { + int isJab; + + VB(("getting cusp info failed\n")); + s->docusp = 0; + + /* ????? Should we use generic cusp information as a fallback ?????? */ + } + + /* Compute source adjustment value */ + for (saval = 0.0, i = 0; i < (sizeof(src_adj)/sizeof(double)-1); i++) + saval += log(src_adj[i]); + saval += src_adj[i]; + + /* For source and dest */ + for (sd = 0; sd < 2; sd++) { + double ta[3] = { 100.0, 0.0, 0.0 }; + double tc[3] = { 0.0, 0.0, 0.0 }; + + /* Compute rotation to rotate/translate so that */ + /* black -> white becomes 0 -> 100 */ + ta[0] *= saval; /* Make source adjustment */ + icmVecRotMat(s->rot[sd], s->cusps[sd][6], s->cusps[sd][7], ta, tc); + + /* And inverse */ + icmVecRotMat(s->irot[sd], ta, tc, s->cusps[sd][6], s->cusps[sd][7]); + + /* Compute a grey */ + if (s->docusp) { + double aL = 0.0; + /* Compute cusp average L value as grey */ + for (k = 0; k < 6; k++) + aL += s->cusps[sd][k][0]; + aL /= 6.0; +//printf("~1 src/dst %d cusp average L %f\n",sd, aL); + + aL = (aL - s->cusps[sd][7][0])/(s->cusps[sd][6][0] - s->cusps[sd][7][0]); /* Param */ + if (aL < 0.0) + aL = 0.0; + else if (aL > 1.0) + aL = 1.0; +//printf("~1 src/dst %d grey param %f\n",sd,aL); + icmBlend3(s->cusps[sd][8], s->cusps[sd][6], s->cusps[sd][7], aL); + + } else { + icmBlend3(s->cusps[sd][8], s->cusps[sd][6], s->cusps[sd][7], 0.5); + } + + /* For white, black and grey */ + icmMul3By3x4(s->cusp_lab[sd][6], s->rot[sd], s->cusps[sd][6]); + icmMul3By3x4(s->cusp_lab[sd][7], s->rot[sd], s->cusps[sd][7]); + icmMul3By3x4(s->cusp_lab[sd][8], s->rot[sd], s->cusps[sd][8]); + + if (!s->docusp) + continue; /* No cusp information */ + + /* For each cusp */ + for (k = 0; k < 6; k++) { + + /* Black/white normalized value */ + icmMul3By3x4(s->cusp_lab[sd][k], s->rot[sd], s->cusps[sd][k]); + + /* Compute LCh value */ + icmLab2LCh(s->cusp_lch[sd][k], s->cusp_lab[sd][k]); + VB(("cusp[%d][%d] %f %f %f LCh %f %f %ff\n", sd, k, s->cusps[sd][k][0], s->cusps[sd][k][1], s->cusps[sd][k][2], s->cusp_lch[sd][k][0], s->cusp_lch[sd][k][1], s->cusp_lch[sd][k][2])); + } + + /* For each pair of cusps */ + for (k = 0; k < 6; k++) { + int m = k < 5 ? k + 1 : 0; + int n; + + /* Plane of grey & 2 cusp points, so as to be able to decide light/dark cone. */ + if (icmPlaneEqn3(s->cusp_pe[sd][k], s->cusp_lab[sd][8], s->cusp_lab[sd][m], + s->cusp_lab[sd][k])) + error("gamut, init_ce: failed to compute plane equation between cusps\n"); + + VB(("dist to white = %f\n",icmPlaneDist3(s->cusp_pe[sd][k], s->cusp_lab[sd][6]))); + VB(("dist to black = %f\n",icmPlaneDist3(s->cusp_pe[sd][k], s->cusp_lab[sd][7]))); + VB(("dist to grey = %f\n",icmPlaneDist3(s->cusp_pe[sd][k], s->cusp_lab[sd][8]))); + VB(("dist to c0 = %f\n",icmPlaneDist3(s->cusp_pe[sd][k], s->cusp_lab[sd][m]))); + VB(("dist to c1 = %f\n",icmPlaneDist3(s->cusp_pe[sd][k], s->cusp_lab[sd][k]))); + + /* For light and dark, create transformation matrix to (src) */ + /* or from (dst) the Baricentric values. The base is always */ + /* the grey point. */ + for (n = 0; n < 2; n++) { + + /* Create from Baricentric matrix */ + icmCpy3(s->cusp_bc[sd][k][n][0], s->cusp_lab[sd][k]); + icmCpy3(s->cusp_bc[sd][k][n][1], s->cusp_lab[sd][m]); + icmCpy3(s->cusp_bc[sd][k][n][2], s->cusp_lab[sd][6 + n]); /* [7] & [8] */ + for (j = 0; j < 3; j++) /* Subtract grey base */ + icmSub3(s->cusp_bc[sd][k][n][j], s->cusp_bc[sd][k][n][j], s->cusp_lab[sd][8]); + + /* Compute matrix transform */ + icmTranspose3x3(s->cusp_bc[sd][k][n], s->cusp_bc[sd][k][n]); + + if (sd == 0) { /* If src, invert matrix */ + if (icmInverse3x3(s->cusp_bc[sd][k][n], s->cusp_bc[sd][k][n]) != 0) + error("gamut, init_ce: failed to invert baricentric matrix\n"); + } + } + } + } + +#ifdef NEVER /* Sanity check */ + + for (k = 0; k < 6; k++) { + double tt[3]; + + comp_ce(s, tt, s->cusps[0][k], NULL); + + VB(("cusp %d, %f %f %f -> %f %f %f, de %f\n", k, cusps[0][k][0], cusps[0][k][1], cusps[0][k][2], tt[0], tt[1], tt[2], icmNorm33(tt, cusps[1][k]))); + + } +#endif /* NEVER */ + +#ifdef NEVER /* Sanity check */ + + { + for (k = 0; k < 9; k++) { + double tt; + + tt = comp_lvc(s, s->cusps[0][k]); + + printf("cusp %d, %f %f %f -> %f\n\n", k, s->cusps[0][k][0], s->cusps[0][k][1], s->cusps[0][k][2], tt); + + } + + /* For light and dark */ + for (sd = 0; sd < 2; sd++) { + /* for each segment */ + for (k = 0; k < 6; k++) { + int m = k < 5 ? k + 1 : 0; + double pos[3], tt; + + pos[0] = pos[1] = pos[2] = 0.0; + + icmAdd3(pos, pos, s->cusps[0][k]); + icmAdd3(pos, pos, s->cusps[0][m]); + icmAdd3(pos, pos, s->cusps[0][6 + sd]); + icmAdd3(pos, pos, s->cusps[0][8]); + icmScale3(pos, pos, 1.0/4.0); + + tt = comp_lvc(s, pos); + + printf("cusps %d & %d, grey %d, %f %f %f -> %f\n\n", k, m, sd, pos[0], pos[1], pos[2], tt); + } + } + } +#endif /* NEVER */ +} + +/* Compute cusp mapping value */ +static void comp_ce( +smthopt *s, /* Context for cusp mapping */ +double out[3], +double in[3], +gammapweights *wt /* If NULL, assume 100% */ +) { + double cw_l = 1.0; + double cw_c = 1.0; + double cw_h = 1.0; + double ccx = 1.0; + + out[0] = in[0]; + out[1] = in[1]; + out[2] = in[2]; + + if (wt != NULL) { + cw_l = wt->c.w.l; + cw_c = wt->c.w.c; + cw_h = wt->c.w.h; + ccx = wt->c.cx; + } + + /* Compute source changes due to any cusp mapping */ + if (s->docusp && (cw_l > 0.0 || cw_c > 0.0 || cw_h > 0.0)) { + double lab[3], lch[3]; /* Normalized source values */ + double bb[3]; /* Baricentric coords */ + double olch[3]; /* Destination transformed LCh source value */ + double mlab[3], mlch[3]; /* Fully mapped value */ + int c0, c1; /* Cusp indexes */ + int ld; /* light/dark index */ + +//printf("\n~1 in = %f %f %f, ccx = %f\n",in[0],in[1],in[2],ccx); + + /* Compute src cusp normalized LCh */ + icmMul3By3x4(lab, s->rot[0], in); + icmLab2LCh(lch, lab); +//printf("~1 lab = %f %f %f, lch = %f %f %f\n",lab[0],lab[1],lab[2],lch[0],lch[1],lch[2]); + + /* Locate the source cusps that this point lies between */ + for (c0 = 0; c0 < 6; c0++) { + double sh, h0, h1; + + sh = lch[2]; + c1 = c0 < 5 ? c0 + 1 : 0; + + h0 = s->cusp_lch[0][c0][2]; + h1 = s->cusp_lch[0][c1][2]; + + if (h1 < h0) { + if (sh < h1) + sh += 360.0; + h1 += 360.0; + } + if (sh >= (h0 - 1e-12) && sh < (h1 + 1e-12)) + break; + } + if (c0 >= 6) /* Assert */ + error("gamut, comp_ce: unable to locate hue %f cusps\n",lch[2]); + + /* See whether this is light or dark */ + ld = icmPlaneDist3(s->cusp_pe[0][c0], lab) >= 0 ? 0 : 1; +//printf("~1 cusp %d, ld %d (dist %f)\n",c0,ld,icmPlaneDist3(s->cusp_pe[0][c0], lab)); + + /* Compute baricentric for input point in simplex */ + icmSub3(bb, lab, s->cusp_lab[0][8]); + icmMulBy3x3(bb, s->cusp_bc[0][c0][ld], bb); + +//printf("~1 bb %f %f %f\n",bb[0],bb[1],bb[2]); + + /* Then compute value for output from baricentric */ + icmMulBy3x3(mlab, s->cusp_bc[1][c0][ld], bb); + icmAdd3(mlab, mlab, s->cusp_lab[1][8]); + icmLab2LCh(mlch, mlab); + +//printf("~1 fully cusp mapped point %f %f %f\n", mlab[0], mlab[1], mlab[2]); + + /* Compute the unchanged source in dest space */ + icmMul3By3x4(olch, s->rot[1], in); + icmLab2LCh(olch, olch); + + /* Then compute weighted output */ + mlch[0] = cw_l * mlch[0] + (1.0 - cw_l) * olch[0]; + mlch[1] = cw_c * mlch[1] + (1.0 - cw_c) * olch[1]; + mlch[1] *= ccx; /* Chroma expansion */ + + if (lch[2] > mlch[2] && (lch[2] - mlch[2]) > 180.0) + mlch[2] += 360.0; + else if (mlch[2] > lch[2] && (mlch[2] - lch[2]) > 180.0) + lch[2] += 360.0; + mlch[2] = cw_c * mlch[2] + (1.0 - cw_c) * lch[2]; + if (mlch[2] >= 360.0) + mlch[2] -= 360.0; + +//printf("~1 weighted cusp mapped point %f %f %f\n", mlch[0], mlch[1], mlch[2]); + icmLCh2Lab(mlch, mlch); + icmMul3By3x4(out, s->irot[1], mlch); +//printf("~1 returning %f %f %f\n", out[0], out[1], out[2]); + } +} + +/* Return a blend factor that measures how close to the white or */ +/* black point the location is. Return 1.0 if as far from the */ +/* point as is grey, 0.0 when at the white or black points. */ +static double comp_naxbf( +smthopt *s, /* Context for cusp mapping */ +double in[3] /* Non-cusp mapped source value */ +) { + double rin[3]; /* Rotated/scaled to neutral axis 0 to 100 */ + double ll; + +//printf("~1 comp_naxbf, %d: in = %f %f %f\n",s->ix, in[0],in[1],in[2]); + + /* Convert to neutral axis 0 to 100 */ + icmMul3By3x4(rin, s->rot[0], in); + +//printf("~1 rotate L %f, white %f grey %f black %f\n",rin[0],s->cusp_lab[0][6][0],s->cusp_lab[0][8][0],s->cusp_lab[0][7][0]); + + if (rin[0] >= s->cusp_lab[0][8][0]) { /* Closer to white */ + ll = icmNorm33(s->cusp_lab[0][6], rin); /* Distance to white */ + ll = 1.0 - ll/(100.0 - s->cusp_lab[0][8][0]); /* Normalized to grey distance */ + } else { /* Closer to black */ + ll = icmNorm33(s->cusp_lab[0][7], rin); /* Distance to black */ + ll = 1.0 - ll/s->cusp_lab[0][8][0]; /* Normalized to grey distance */ + } + if (ll < 0.0) + ll = 0.0; + else if (ll > 1.0) + ll = 1.0; + + /* Weight so that it goes to 0.0 close to W & B */ + ll = sqrt(1.0 - ll); + +//printf("~1 returning ll %f\n",ll); + return ll; +} + +/* Return a value suitable for blending between the wl, gl and bl L dominance values */ +/* The value is a linear blend value, 0.0 at cusp local grey, 1.0 at white L value */ +/* and -1.0 at black L value. */ +static double comp_lvc( +smthopt *s, /* Context for cusp mapping */ +double in[3] /* Non-cusp mapped source value */ +) { + double Lg; + double ll; + +//printf("~1 comp_lvc, %d: in = %f %f %f\n",s->ix, in[0],in[1],in[2]); + + /* Compute the cusp local grey value. */ + if (s->docusp) { + double lab[3], lch[3]; /* Normalized source values */ + double bb[3]; /* Baricentric coords */ + int c0, c1; /* Cusp indexes */ + int ld; /* light/dark index */ + + /* Compute src cusp normalized LCh */ + icmMul3By3x4(lab, s->rot[0], in); + icmLab2LCh(lch, lab); +//printf("~1 lab = %f %f %f, lch = %f %f %f\n",lab[0],lab[1],lab[2],lch[0],lch[1],lch[2]); + + /* Locate the source cusps that this point lies between */ + for (c0 = 0; c0 < 6; c0++) { + double sh, h0, h1; + + sh = lch[2]; + c1 = c0 < 5 ? c0 + 1 : 0; + + h0 = s->cusp_lch[0][c0][2]; + h1 = s->cusp_lch[0][c1][2]; + + if (h1 < h0) { + if (sh < h1) + sh += 360.0; + h1 += 360.0; + } + if (sh >= (h0 - 1e-12) && sh < (h1 + 1e-12)) + break; + } + if (c0 >= 6) /* Assert */ + error("gamut, comp_lvc: unable to locate hue %f cusps\n",lch[2]); + + /* See whether this is light or dark */ + ld = icmPlaneDist3(s->cusp_pe[0][c0], lab) >= 0 ? 0 : 1; +//printf("~1 cusp %d, ld %d (dist %f)\n",c0,ld,icmPlaneDist3(s->cusp_pe[0][c0], lab)); + + /* Compute baricentric for input point in simplex */ + icmSub3(bb, lab, s->cusp_lab[0][8]); + icmMulBy3x3(bb, s->cusp_bc[0][c0][ld], bb); + +//printf("~1 baricentric %f %f %f\n",bb[0],bb[1],bb[2]); + + /* Compute the grey level */ + Lg = s->cusps[0][8][0] + + bb[0] * (s->cusps[0][c0][0] - s->cusps[0][8][0]) + + bb[1] * (s->cusps[0][c1][0] - s->cusps[0][8][0]); + + } else { + /* Non-cusp sensitive grey L */ + Lg = s->cusps[0][8][0]; + } +//printf("~1 grey = %f\n",Lg); + + if (in[0] > Lg) { + ll = (in[0] - Lg)/(s->cusps[0][6][0] - Lg); + + } else { + ll = -(in[0] - Lg)/(s->cusps[0][7][0] - Lg); + } +//printf("~1 returnin ll %f\n",ll); + return ll; +} + +static double invfunc( +void *fdata, +double *tp +) { + smthopt *s = (smthopt *)fdata; + double cv[3]; /* Converted value */ + double tt, rv = 0.0; + + comp_ce(s, cv, tp, s->wt); + + tt = s->tv[0] - cv[0]; + rv += tt * tt; + tt = s->tv[1] - cv[1]; + rv += tt * tt; + tt = s->tv[2] - cv[2]; + rv += tt * tt; + +//printf("~1 rv %f from %f %f %f -> %f %f %f\n",rv,tp[0],tp[1],tp[2],cv[0],cv[1],cv[2]); + return rv; +} + +/* Inverse of com_ce. We do this by inverting com_ce numerically (slow) */ +static void inv_comp_ce( +smthopt *s, /* Context for cusp mapping */ +double out[3], +double in[3], +gammapweights *wt /* If NULL, assume 100% */ +) { + double ss[3] = { 20.0, 20.0, 20.0 }; /* search area */ + double tp[3], rv; + + s->tv[0] = tp[0] = in[0]; + s->tv[1] = tp[1] = in[1]; + s->tv[2] = tp[2] = in[2]; + s->wt = wt; + + /* Optimise the point */ + if (powell(&rv, 3, tp, ss, 0.001, 1000, invfunc, (void *)s, NULL, NULL) != 0) { + error("gammap::nearsmth: inv_comp_ce powell failed on %f %f %f\n", in[0], in[1], in[2]); + } + +//printf("~1 inv_comp_ce: %f %f %f -> %f %f %f\n", s->tv[0], s->tv[1], s->tv[2], tp[0], tp[1], tp[2]); +//comp_ce(s, out, tp, wt); +//printf("~1 check: %f %f %f, DE %f\n", out[0], out[1], out[2], icmNorm33(s->tv,out)); + + out[0] = tp[0]; + out[1] = tp[1]; + out[2] = tp[2]; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* A set of functions to help handle the weighting configuration */ + +/* Copy non-negative values from one set of weights to another */ +void near_wcopy( +gammapweights *dst, +gammapweights *src +) { + +#define NSCOPY(xxx) dst->xxx = src->xxx >= 0.0 ? src->xxx : dst->xxx +//#define NSCOPY(xxx) if (src->xxx >= 0.0) { \ +// printf("Setting %s to %f\n",#xxx, src->xxx); \ +// dst->xxx = src->xxx; \ +// } + + NSCOPY(c.w.l); + NSCOPY(c.w.c); + NSCOPY(c.w.h); + NSCOPY(c.cx); + + NSCOPY(l.o); + NSCOPY(l.h); + NSCOPY(l.l); + + NSCOPY(a.o); + NSCOPY(a.h); + NSCOPY(a.wl); + NSCOPY(a.gl); + NSCOPY(a.bl); + + NSCOPY(a.wlth); + NSCOPY(a.blpow); + + NSCOPY(a.lxpow); + NSCOPY(a.lxthr); + + NSCOPY(r.rdl); + NSCOPY(r.rdh); + + NSCOPY(d.co); + NSCOPY(d.xo); + + NSCOPY(f.x); +#undef NSCOPY +} + +/* Blend a two groups of individual weights into one, given two weightings */ +void near_wblend( +gammapweights *dst, +gammapweights *src1, double wgt1, +gammapweights *src2, double wgt2 +) { + +#define NSBLEND(xxx) dst->xxx = wgt1 * src1->xxx + wgt2 * src2->xxx + + NSBLEND(c.w.l); + NSBLEND(c.w.c); + NSBLEND(c.w.h); + NSBLEND(c.cx); + + NSBLEND(l.o); + NSBLEND(l.h); + NSBLEND(l.l); + + NSBLEND(a.o); + NSBLEND(a.h); + NSBLEND(a.wl); + NSBLEND(a.gl); + NSBLEND(a.bl); + + NSBLEND(a.wlth); + NSBLEND(a.blpow); + + NSBLEND(a.lxpow); + NSBLEND(a.lxthr); + + NSBLEND(r.rdl); + NSBLEND(r.rdh); + + NSBLEND(d.co); + NSBLEND(d.xo); + + NSBLEND(f.x); +#undef NSBLEND +} + +/* Expand the compact form of weights into the explicit form. */ +/* The explicit form is light and dark of red, yellow, green, cyan, blue, magenta & neutral*/ +/* Return nz on error */ +int expand_weights(gammapweights out[14], gammapweights *in) { + int i, j; + + /* Set the usage of each slot */ + out[0].ch = gmm_light_red; + out[1].ch = gmm_light_yellow; + out[2].ch = gmm_light_green; + out[3].ch = gmm_light_cyan; + out[4].ch = gmm_light_blue; + out[5].ch = gmm_light_magenta; + out[6].ch = gmm_light_neutral; + + out[7].ch = gmm_dark_red; + out[8].ch = gmm_dark_yellow; + out[9].ch = gmm_dark_green; + out[10].ch = gmm_dark_cyan; + out[11].ch = gmm_dark_blue; + out[12].ch = gmm_dark_magenta; + out[13].ch = gmm_dark_neutral; + +//printf("\n~1 expand weights called\n"); + + /* mark output so we can recognise having been set or not */ + for (i = 0; i < 14; i++) + out[i].set = 0; + + /* Expand the compact form to explicit. */ + + /* First is default */ + for (i = 0; in[i].ch != gmm_end; i++) { + if (in[i].ch == gmm_end) + break; + if (in[i].ch == gmm_ignore) + continue; + + if (in[i].ch == gmm_default) { + for (j = 0; j < 14; j++) { +//printf("~1 Setting %d 0x%x with 0x%x (default)\n",j,out[j].ch,in[i].ch); + if ((in[i].ch & out[j].ch) == out[j].ch) { + near_wcopy(&out[j], &in[i]); + out[j].set = 1; + } + } + } + } + + /* Then light or dark */ + for (i = 0; in[i].ch != gmm_end; i++) { + if (in[i].ch == gmm_end) + break; + if (in[i].ch == gmm_ignore) + continue; + + if (in[i].ch == gmm_light_colors + || in[i].ch == gmm_dark_colors) { + for (j = 0; j < 14; j++) { + if ((in[i].ch & out[j].ch) == out[j].ch) { +//printf("~1 Setting %d 0x%x with 0x%x (light or dark)\n",j,out[j].ch,in[i].ch); + near_wcopy(&out[j], &in[i]); + out[j].set = 1; + } + } + } + } + + /* Then light and dark colors */ + for (i = 0; in[i].ch != gmm_end; i++) { + if (in[i].ch == gmm_end) + break; + if (in[i].ch == gmm_ignore) + continue; + + if ((in[i].ch & gmc_l_d) == gmc_l_d + && (in[i].ch & gmc_colors) != gmc_colors) { + for (j = 0; j < 14; j++) { + if ((in[i].ch & out[j].ch) == out[j].ch) { +//printf("~1 Setting %d 0x%x with 0x%x (light and dark color)\n",j,out[j].ch,in[i].ch); + near_wcopy(&out[j], &in[i]); + out[j].set = 1; + } + } + } + } + + /* Last pass is light or dark colors */ + for (i = 0; in[i].ch != gmm_end; i++) { + if (in[i].ch == gmm_end) + break; + if (in[i].ch == gmm_ignore) + continue; + + if (((in[i].ch & gmc_l_d) == gmc_light + || (in[i].ch & gmc_l_d) == gmc_dark) + && (in[i].ch & gmc_colors) != gmc_colors) { + for (j = 0; j < 14; j++) { + if ((in[i].ch & out[j].ch) == out[j].ch) { +//printf("~1 Setting %d 0x%x with 0x%x (light or dark color)\n",j,out[j].ch,in[i].ch); + near_wcopy(&out[j], &in[i]); + out[j].set = 1; + } + } + } + } + + /* Check every slot has been set */ + for (i = 0; i < 14; i++) { + if (out[i].set == 0) { +//printf("~1 set %d hasn't been initialized\n",i); + return 1; + } + } + return 0; +} + +/* Tweak weights acording to extra cmy cusp flags or rel override */ +void tweak_weights(gammapweights out[14], int dst_cmymap, int rel_oride) { + int i; + + for (i = 0; i < 14; i++) { + if (((dst_cmymap & 0x1) && (out[i].ch & gmc_cyan)) + || ((dst_cmymap & 0x2) && (out[i].ch & gmc_magenta)) + || ((dst_cmymap & 0x4) && (out[i].ch & gmc_yellow))) { +//printf("~1 Setting %d 0x%x to 100% cusp map\n",i,out[i].ch); + out[i].c.w.l = 1.0; /* 100% mapping */ + out[i].c.w.c = 1.0; + out[i].c.w.h = 1.0; + out[i].c.cx = 1.0; /* No expansion */ + } + + if (rel_oride == 1) { /* A high saturation "clip" like mapping */ + out[i].r.rdl = 1.0; /* No relative neighbourhood */ + out[i].r.rdh = 1.0; /* No relative neighbourhood */ + out[i].d.co = 0.0; /* No depth weighting */ + out[i].d.xo = 0.0; /* No depth weighting */ + + } else if (rel_oride == 2) { /* A maximal feature preserving mapping */ + out[i].r.rdl *= 1.6; /* Extra neighbourhood size */ + out[i].r.rdh *= 1.6; /* Extra neighbourhood size */ + } + } +} + +/* Blend a two expanded groups of individual weights into one */ +void near_xwblend( +gammapweights *dst, +gammapweights *src1, double wgt1, +gammapweights *src2, double wgt2 +) { + int i; + for (i = 0; i < 14; i++) + near_wblend(&dst[i], &src1[i], wgt1, &src2[i], wgt2); +} + +/* Convert overall, hue dom & l dom to iweight */ +static void comp_iweight(iweight *iw, double o, double h, double l) { + double c, lc; + + if (h < 0.0) + h = 0.0; + else if (h > 1.0) + h = 1.0; + + if (l < 0.0) + l = 0.0; + else if (l > 1.0) + l = 1.0; + + lc = 1.0 - h; + c = (1.0 - l) * lc; + l = l * lc; + + o /= sqrt(l * l + c * c + h * h); + + iw->l = o * l; + iw->c = o * c; + iw->h = o * h; +} + +/* Given a point location, return the interpolated weighting values at that point. */ +/* (Typically non-cusp mapped source location assumed, and source gamut cusps used) */ +/* (Assume init_ce() has been called to setip smthopt!) */ +void interp_xweights(gamut *gam, gammapweights *out, double pos[3], + gammapweights in[14], smthopt *s, int cvec) { + double h, JCh[3], tmp[3]; + int li, ui; /* The two hue indexes the color is between */ + double lh, uh; /* Lower/upper hue of two colors */ + double lw, uw; /* Lower/upper blend values */ + double cusps[6][3]; + gammapweights light, dark; + + /* Convert to polar */ + icmLab2LCh(JCh, pos); + + if (gam->getcusps(gam, cusps) != 0) { /* Failed */ + int isJab = gam->isJab ? 1 : 0; + + /* Figure out what hextant we're between using generic cusps */ + for (li = 0; li < 6; li++) { + ui = li < 5 ? li + 1 : 0; + h = JCh[2]; + + lh = gam_hues[isJab][li]; /* use generic ones */ + uh = gam_hues[isJab][ui]; + if (uh < lh) { + if (h < uh) + h += 360.0; + uh += 360.0; + } + if (h >= (lh - 1e-12) && h < (uh + 1e-12)) + break; + } + + } else { + + /* Locate the source cusps that this point lies between */ + for (li = 0; li < 6; li++) { + double tt[3]; + ui = li < 5 ? li + 1 : 0; + h = JCh[2]; + + icmLab2LCh(tt, cusps[li]); + lh = tt[2]; + icmLab2LCh(tt, cusps[ui]); + uh = tt[2]; + + if (uh < lh) { + if (h < uh) + h += 360.0; + uh += 360.0; + } + if (h >= (lh - 1e-12) && h < (uh + 1e-12)) + break; + } + } + if (li >= 6) /* Assert */ + error("gamut, interp_xweights: unable to locate hue %f cusps\n",JCh[2]); + + /* Compute hue angle blend weights */ + uw = (h - lh)/(uh - lh); + if (uw < 0.0) + uw = 0.0; + else if (uw > 1.0) + uw = 1.0; + uw = uw * uw * (3.0 - 2.0 * uw); /* Apply spline to smooth interpolation */ + lw = (1.0 - uw); + + /* Blend weights at the two hues */ + near_wblend(&light, &in[li], lw, &in[ui], uw); + near_wblend(&dark, &in[7 + li], lw, &in[7 + ui], uw); + + /* If we're close to the center, blend to the neutral weight */ + if (JCh[1] < NEUTRAL_C) { + lw = (NEUTRAL_C - JCh[1])/NEUTRAL_C; + uw = (1.0 - lw); + near_wblend(&light, &in[6], lw, &light, uw); + near_wblend(&dark, &in[7 + 6], lw, &dark, uw); + } + + /* Figure out where we are between light and dark, */ + /* and create blend between their weightings */ + uw = (JCh[0] - DARK_L)/(LIGHT_L - DARK_L); + if (uw > 1.0) + uw = 1.0; + else if (uw < 0.0) + uw = 0.0; + uw = uw * uw * (3.0 - 2.0 * uw); /* Apply spline to smooth interpolation */ + lw = (1.0 - uw); + near_wblend(out, &dark, lw, &light, uw); + + /* Convert radial dominance weights into raw weights */ + comp_iweight(&out->rl, out->l.o, out->l.h, out->l.l); + +//printf("~1 %d: src %f %f %f (cvec %d)\n",s->ix, pos[0],pos[1],pos[2],cvec); + /* Compute l dominance value vs. closness to white or black point */ + { + double wl, gl, bl, uw, l; + + /* Closness to white and black points */ + uw = comp_lvc(s, pos); + +//printf("~1 uw = %f\n",uw); + if (uw >= 0) { + + /* Scale to threshold */ + if (uw > (1.0 - out->a.wlth)) + uw = (uw - 1.0 + out->a.wlth)/out->a.wlth; + else + uw = 0.0; +//printf("~1 white, thresholded uw %f\n",uw); + + /* Blend in log ratio space */ + wl = log((1.0 - out->a.wl + 1e-5)/(out->a.wl + 1e-5)); + gl = log((1.0 - out->a.gl + 1e-5)/(out->a.gl + 1e-5)); + l = exp(uw * wl + (1.0 - uw) * gl); + l = ((1.0 - l) * 1e-5 + 1.0)/(l + 1.0); + + } else { + uw = -uw; + + /* Apply power */ + uw = pow(uw, out->a.blpow); +//printf("~1 black with power uw %f\n",uw); + + /* Blend in log ratio space */ + gl = log((1.0 - out->a.gl + 1e-5)/(out->a.gl + 1e-5)); + bl = log((1.0 - out->a.bl + 1e-5)/(out->a.bl + 1e-5)); + l = exp(uw * bl + (1.0 - uw) * gl); + l = ((1.0 - l) * 1e-5 + 1.0)/(l + 1.0); + } +//printf("~1 wl %f, gl %f, bl %f -> %f\n",out->a.wl,out->a.gl,out->a.bl,l); + + /* Convert absolute dominance weights into raw weights */ + comp_iweight(&out->ra, out->a.o, out->a.h, l); +//printf("~1 l %f, h %f, ra l %f c %f h %f\n\n", l, out->a.h, out->ra.l, out->ra.c, out->ra.h); + } +} + +/* Callback used by compdstgamut() to establish the expected compression */ +/* mapping direction. */ +static void cvect( +void *cntx, /* smthopt * */ +double *p2, /* Return point displaced from p1 in desired direction */ +double *p1 /* Given point */ +) { + double vv, gv[3], lv[3]; + smthopt *s = (smthopt *)cntx; + gammapweights out; + + interp_xweights(s->sgam, &out, p1, s->xwh, s, 1); + +//printf("~1 at %f %f %f, lch weight %f %f %f\n", p1[0], p1[1], p1[2], out.ra.l, out.ra.c, out.ra.h); + /* Now we need to convert the absolute weighting out.ra into a vector */ + /* We do this in a very simple minded fashion. The hue weighting is ignored, */ + /* because we assume a direction towards the neutral axis. The C weight is */ + /* assumed to be the weight towards the grey point, while the L weight */ + /* assumed to be the weight towards the point on the neutral axis with */ + /* the same L value. */ + + /* Parameter along neutral axis black to white */ + vv = (p1[0] - s->cusps[0][7][0])/(s->cusps[0][6][0] - s->cusps[0][7][0]); + /* lv is point at same L on neutral axis */ + lv[0] = p1[0]; + lv[1] = vv * (s->cusps[0][6][1] - s->cusps[0][7][1]) + s->cusps[0][7][1]; + lv[2] = vv * (s->cusps[0][6][2] - s->cusps[0][7][2]) + s->cusps[0][7][2]; + icmSub3(lv, lv, p1); /* Vector */ + icmNormalize3(lv, lv, out.ra.l); + + icmSub3(gv, s->cusps[0][8], p1); /* Grey vector */ + icmNormalize3(gv, gv, out.ra.c); + + icmAdd3(p2, gv, p1); + icmAdd3(p2, lv, p2); /* Combined destination */ +//printf("~1 p2 %f %f %f\n", p2[0], p2[1], p2[2]); +} + +/* Shrink function */ +static void doshrink(smthopt *s, double *out, double *in, double shrink) { + double rad, len, p2[3]; + + cvect((void *)s, p2, in); /* Get shrink direction */ + + /* Conservative radius of point */ + rad = sqrt(in[1] * in[1] + in[2] * in[2]); + + len = shrink; + if (rad < (2.0 * shrink)) + len = rad * 0.5; + + icmNormalize33(out, p2, in, len); +} + +/* Convenience function. Given a mapping vector, return the */ +/* intersection with the given gamut that is in the mapping */ +/* direction. Return NZ if no intersection */ +static int vintersect( +gamut *g, +int *p1out, /* Return nz if p1 is outside the gamut */ +double isec[3], /* Return intersection point */ +double p1[3], /* First point */ +double p2[3] /* Second point */ +) { + gispnt lp[40]; + int ll, i, bi; + + if ((ll = g->vector_isectns(g, p1, p2, lp, 40)) == 0) + return 1; + + /* Locate the segment or non-segment the source lies in */ + for (bi = -1, i = 0; i < ll; i += 2) { + + if ((i == 0 || lp[i-1].pv < 0.0) /* p1 is outside gamut */ + && lp[i].pv >= -1e-2) { + bi = i; + if (p1out != NULL) + *p1out = 1; + break; + } + + if (lp[i].pv <= 0.0 /* p1 is inside gamut */ + && lp[i+1].pv >= -1e-2) { + bi = i+1; + if (p1out != NULL) + *p1out = 0; + break; + } + } + if (bi < 0) + return 1; + + if (isec != NULL) + icmCpy3(isec, lp[bi].ip); + + return 0; +} + +/* Convenience function. Given a point and an inwards mapping vector, */ +/* if the point is within the gamut, return the first intersection in */ +/* the opposite to vector direction. If the point is outside the gamut, */ +/* return the first intersction in the vector direction. */ +/* Return NZ if no intersection */ +static int vintersect2( +gamut *g, +int *p1out, /* Return nz if p1 is outside the gamut */ +double isec[3], /* Return intersection point */ +double vec[3], /* Vector */ +double p1[3] /* Point */ +) { + gispnt lp[40]; + double p2[3]; + int ll, i, bi; + + icmAdd3(p2, p1, vec); + + if ((ll = g->vector_isectns(g, p1, p2, lp, 40)) == 0) + return 1; + + /* Locate the segment or non-segment the source lies in */ + for (bi = -1, i = 0; i < ll; i += 2) { + + if ((i == 0 || lp[i-1].pv < 0.0) /* p1 is outside gamut, */ + && lp[i].pv >= 0.0) { /* so look in +ve pv direction. */ + bi = i; + if (p1out != NULL) + *p1out = 1; + break; + } + + if (lp[i].pv <= 0.0 /* p1 is inside gamut, */ + && lp[i+1].pv >= 0.0) { /* so look in -ve pv direction. */ + bi = i; + if (p1out != NULL) + *p1out = 0; + break; + } + } + if (bi < 0) + return 1; + + if (isec != NULL) + icmCpy3(isec, lp[bi].ip); + + return 0; +} + + +/* ============================================ */ + +/* Return the maximum number of points that will be generated */ +int near_smooth_np( + gamut *sc_gam, /* Source colorspace gamut */ + gamut *si_gam, /* Source image gamut (== sc_gam if none) */ + gamut *d_gam, /* Destination colorspace gamut */ + double xvra /* Extra vertex ratio */ +) { + gamut *p_gam; /* Gamut used for points - either source colorspace or image */ + int ntpts, nmpts, nspts, nipts, ndpts; + + nspts = sc_gam->nverts(sc_gam); + nipts = si_gam->nverts(si_gam); + ndpts = d_gam->nverts(d_gam); + p_gam = sc_gam; + + /* Target number of points is max of any gamut */ + ntpts = nspts > nipts ? nspts : nipts; + ntpts = ntpts > ndpts ? ntpts : ndpts; + ntpts = (int)(ntpts * xvra + 0.5); + + /* Use image gamut if it exists */ + if (nspts < nipts || si_gam != sc_gam) { + nspts = nipts; /* Use image gamut instead */ + p_gam = si_gam; + } + xvra = ntpts/(double)nspts; + nmpts = p_gam->nssverts(p_gam, xvra); /* Stratified Sampling source points */ + + return nmpts; +} + + +/* ============================================ */ +/* Return a list of points. Free list after use */ +/* Return NULL on error */ +nearsmth *near_smooth( +int verb, /* Verbose flag */ +int *npp, /* Return the actual number of points returned */ +gamut *sc_gam, /* Source colorspace gamut - uses cusp info if availablle */ +gamut *si_gam, /* Source image gamut (== sc_gam if none), just used for surface. */ +gamut *d_gam, /* Destination colorspace gamut */ +int src_kbp, /* Use K only black point as src gamut black point */ +int dst_kbp, /* Use K only black point as dst gamut black point */ +double d_bp[3], /* Override destination target black point - may be NULL */ +gammapweights xwh[14],/* Structure holding expanded hextant weightings */ +double gamcknf, /* Gamut compression knee factor, 0.0 - 1.0 */ +double gamxknf, /* Gamut expansion knee factor, 0.0 - 1.0 */ +int usecomp, /* Flag indicating whether smoothed compressed value will be used */ +int useexp, /* Flag indicating whether smoothed expanded value will be used */ +double xvra, /* Extra number of vertexes ratio */ +int mapres, /* Grid res for 3D RSPL */ +double mapsmooth, /* Target smoothing for 3D RSPL */ +datai map_il, /* Preliminary rspl input range */ +datai map_ih, +datao map_ol, /* Preliminary rspl output range */ +datao map_oh +) { + smthopt opts; /* optimisation and cusp mapping context */ + int ix, i, j, k; + gamut *p_gam; /* Gamut used for points - either source colorspace or image */ + gamut *sci_gam; /* Intersection of src and img gamut gamut */ + gamut *di_gam; /* Modified destination gamut suitable for mapping from sci_gam. */ + /* If just compression, this is the smaller of sci_gam and d_gam. */ + /* If just expansion, this is the sci_gam expanded by d_gam - sc_gam. */ + /* If both exp & comp, then where + d_gam is outside sci_gam this is sci_gam expanded by d_gam - sc_gam + else it is the smaller of sci_gam and d_gam */ + gamut *nedi_gam;/* Same as above, but not expanded. */ + int nmpts; /* Number of mapping gamut points */ + nearsmth *smp; /* Absolute delta E weighting */ + int pass; + int it; + rspl *evectmap = NULL; /* evector map */ + double codf; /* Itteration overshoot/damping factor */ + double mxmv; /* Maximum a point gets moved */ + int nmxmv; /* Number of maxmoves less than stopping threshold */ + + /* Check gamuts are compatible */ + if (sc_gam->compatible(sc_gam, d_gam) == 0 + || (si_gam != NULL && sc_gam->compatible(sc_gam, si_gam) == 0)) { + fprintf(stderr,"gamut map: Gamuts aren't compatible\n"); + *npp = 0; + return NULL; + } + + { + int ntpts, nspts, nipts, ndpts; + + nspts = sc_gam->nverts(sc_gam); + nipts = si_gam->nverts(si_gam); + ndpts = d_gam->nverts(d_gam); + p_gam = sc_gam; + + /* Target number of points is max of any gamut */ + ntpts = nspts > nipts ? nspts : nipts; + ntpts = ntpts > ndpts ? ntpts : ndpts; + ntpts = (int)(ntpts * xvra + 0.5); + + /* Use image gamut if it exists */ + if (nspts < nipts || si_gam != sc_gam) { + nspts = nipts; /* Use image gamut instead */ + p_gam = si_gam; + } + xvra = ntpts/(double)nspts; + nmpts = p_gam->nssverts(p_gam, xvra); /* Stratified Sampling source points */ + + if (verb) printf("Vertex count: mult. = %f, src %d, img %d dst %d, target %d\n", + xvra,nspts,nipts,ndpts,nmpts); + } + + /* Setup opts structure */ + opts.useexp = useexp; /* Expansion used ? */ + opts.debug = 0; /* No debug powell() failure */ + opts.xwh = xwh; /* Weightings */ + opts.sgam = sc_gam; /* Source colorspace gamut */ + + /* Setup source & dest neutral axis transform if white/black available. */ + /* If cusps are available, also figure out the transformations */ + /* needed to map source cusps to destination cusps */ + init_ce(&opts, sc_gam, d_gam, src_kbp, dst_kbp, d_bp); + + /* Allocate our guide points */ + if ((smp = (nearsmth *)calloc(nmpts, sizeof(nearsmth))) == NULL) { + fprintf(stderr,"gamut map: Malloc of near smooth points failed\n"); + *npp = 0; + return NULL; + } + + /* Create a source gamut surface that is the intersection of the src colorspace */ + /* and image gamut, in case (for some strange reason) the image gamut. */ + /* exceeds the source colorspace size. */ + sci_gam = sc_gam; /* Alias to source space gamut */ + if (si_gam != sc_gam) { + if ((sci_gam = new_gamut(0.0, 0, 0)) == NULL) { + fprintf(stderr,"gamut map: new_gamut failed\n"); + *npp = 0; + return NULL; + } + sci_gam->intersect(sci_gam, sc_gam, si_gam); +#ifdef SAVE_VRMLS + printf("###### gamut/nearsmth.c: writing diagnostic sci_gam.wrl and di_gam.wrl\n"); + sci_gam->write_vrml(sci_gam, "sci_gam.wrl", 1, 0); +#endif + } + + di_gam = sci_gam; /* Default no compress or expand */ + if (usecomp || useexp) { + if ((di_gam = new_gamut(0.0, 0, 0)) == NULL) { + fprintf(stderr,"gamut map: new_gamut failed\n"); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + *npp = 0; + return NULL; + } + if (useexp) { + /* Non-expand di_gam for testing double back img points against */ + if ((nedi_gam = new_gamut(0.0, 0, 0)) == NULL) { + fprintf(stderr,"gamut map: new_gamut failed\n"); + di_gam->del(di_gam); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + *npp = 0; + return NULL; + } + } else { + nedi_gam = di_gam; + } + + /* Initialise this gamut with a destination mapping gamut. */ + /* This will be the smaller of the image and destination gamut if compressing, */ + /* and will be expanded by the (dst - src) if expanding. */ + di_gam->compdstgamut(di_gam, sci_gam, sc_gam, d_gam, usecomp, useexp, nedi_gam, + cvect, &opts); + } + +#ifdef SAVE_VRMLS + di_gam->write_vrml(di_gam, "di_gam.wrl", 1, 0); +#endif + + /* Create a list of the mapping guide points, setup for a null mapping */ + VA(("Creating the mapping guide point list\n")); + for (ix = i = 0; i < nmpts; i++) { + double imv[3], imr; /* Image gamut source point and radius */ + double inorm[3]; /* Normal of image gamut surface at src point */ + + /* Get the source color/image space vertex value we are going */ + /* to use as a sample point. */ + if ((ix = p_gam->getssvert(p_gam, &imr, imv, inorm, ix)) < 0) { + break; + } +//printf("~1 got point %d out of %d\n",i+1,nmpts); + + if (p_gam != sc_gam) { /* If src colorspace point, map to img gamut surface */ + imr = sci_gam->radial(sci_gam, imv, imv); + } + + /* If point is within non-expanded modified destination gamut, */ + /* then it is a "double back" image point, and should be ignored. */ + if (nedi_gam->radial(nedi_gam, NULL, imv) > (imr + 1e-4)) { + VB(("Rejecting point %d because it's inside destination\n",i)); + i--; + continue; + } + + /* Lookup radialy equivalent point on modified destination gamut, */ + /* in case we need it for compression or expansion */ + smp[i].drr = di_gam->radial(di_gam, smp[i].drv, imv); + + /* Default setup a null mapping of source image space point to source image point */ + smp[i].vflag = smp[i].gflag = 0; + smp[i].dr = smp[i].sr = smp[i]._sr = imr; + smp[i].dv[0] = smp[i].sv[0] = smp[i]._sv[0] = imv[0]; + smp[i].dv[1] = smp[i].sv[1] = smp[i]._sv[1] = imv[1]; + smp[i].dv[2] = smp[i].sv[2] = smp[i]._sv[2] = imv[2]; + smp[i].sgam = sci_gam; + smp[i].dgam = sci_gam; + + VB(("In Src %d = %f %f %f\n",i,smp[i].sv[0],smp[i].sv[1],smp[i].sv[2])); + + /* If we're going to comp. or exp., check that the guide vertex is not */ + /* on the wrong side of the image gamut, due to the it being */ + /* a small subset of the source colorspace, displaced to one side. */ + /* Because of the gamut convexity limitations, this amounts */ + /* to the source surface at the vertex being in the direction */ + /* of the center. */ + if (usecomp != 0 || useexp != 0) { + double mv[3], ml; /* Radial inward mapping vector */ + double dir; + + icmSub3(mv, sci_gam->cent, smp[i].sv); /* Vector to center */ + ml = icmNorm3(mv); /* It's length */ + + if (ml > 0.001) { + dir = icmDot3(mv, inorm); /* Compare to normal of src triangle */ +//printf("~1 ix %d, dir = %f, dir/len = %f\n",i,dir, dir/ml); + dir /= ml; + if (dir < 0.02) { /* If very shallow */ +//printf("~1 rejecting point %d because it's oblique\n",i); + VB(("Rejecting point %d because it's oblique\n",i)); + i--; + continue; + } + } + } + + /* Set some default extra guide point values */ + smp[i].anv[0] = smp[i].aodv[0] = smp[i].dv[0]; + smp[i].anv[1] = smp[i].aodv[1] = smp[i].dv[1]; + smp[i].anv[2] = smp[i].aodv[2] = smp[i].dv[2]; + + VB(("Src %d = %f %f %f\n",i,smp[i].sv[0],smp[i].sv[1],smp[i].sv[2])); + VB(("Dst %d = %f %f %f\n",i,smp[i].dv[0],smp[i].dv[1],smp[i].dv[2])); + } + nmpts = i; /* Number of points after rejecting any */ + *npp = nmpts; + + /* Don't need this anymore */ + if (nedi_gam != di_gam) + nedi_gam->del(nedi_gam); + nedi_gam = NULL; + + /* If nothing to be compressed or expanded, then return */ + if (usecomp == 0 && useexp == 0) { + VB(("Neither compression nor expansion defined\n")); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + return smp; + } + + /* Set the parameter weights for each point */ + for (i = 0; i < nmpts; i++) { + opts.ix = i; /* Point in question */ + opts.p = &smp[i]; + + /* Determine the parameter weighting for this point */ + interp_xweights(opts.sgam, &smp[i].wt, smp[i]._sv, opts.xwh, &opts, 0); + } + + VA(("Setting up cusp rotated compression or expansion mappings\n")); + VB(("rimv = Cusp rotated cspace/image gamut source point\n")); + VB(("imv = cspace/image gamut source point\n")); + VB(("drv = Destination space radial point and radius \n")); + + /* Setup the cusp rotated compression or expansion mappings */ + for (i = 0; i < nmpts; i++) { + double imv[3], imr; /* cspace/image gamut source point and radius */ + double rimv[3], rimr; /* Cusp rotated cspace/image gamut source point and radius */ + + opts.ix = i; /* Point in question */ + opts.p = &smp[i]; + + /* Grab the source image point */ + imr = smp[i]._sr; + imv[0] = smp[i]._sv[0]; + imv[1] = smp[i]._sv[1]; + imv[2] = smp[i]._sv[2]; + + /* Compute the cusp rotated version of the cspace/image points */ + /* Note that we're not elevating yet! */ + comp_ce(&opts, rimv, imv, &smp[i].wt); + VB(("%f de, ix %d: cusp mapped %f %f %f -> %f %f %f\n", icmNorm33(rimv,imv), i, imv[0], imv[1], imv[2], rimv[0], rimv[1], rimv[2])); + rimr = icmNorm33(rimv, sci_gam->cent); + + /* Default setup a no compress or expand mapping of */ + /* source space/image point to modified destination gamut. */ + smp[i].sr = rimr; + smp[i].sv[0] = rimv[0]; /* Temporary rotated src point */ + smp[i].sv[1] = rimv[1]; + smp[i].sv[2] = rimv[2]; + smp[i].sgam = sci_gam; + smp[i].dgam = di_gam; + + VB(("\n")); + VB(("point %d:, rimv = %f %f %f, rimr = %f\n",i,rimv[0],rimv[1],rimv[2],rimr)); + VB(("point %d:, imv = %f %f %f, imr = %f\n",i,imv[0],imv[1],imv[2],imr)); + VB(("point %d:, drv = %f %f %f, drr = %f\n",i,smp[i].drv[0],smp[i].drv[1],smp[i].drv[2],smp[i].drr)); + + /* Set a starting point for the optimisation */ + smp[i].dgam->nearest(smp[i].dgam, smp[i].dv, smp[i].sv); + smp[i].dr = icmNorm33(smp[i].dv, smp[i].dgam->cent); + + /* Re-lookup radialy equivalent point on destination gamut, */ + /* to match rotated/elevated source */ + smp[i].drr = smp[i].dgam->radial(smp[i].dgam, smp[i].drv, smp[i].sv); + + /* A default average neighbour value */ + smp[i].anv[0] = smp[i].drv[0]; + smp[i].anv[1] = smp[i].drv[1]; + smp[i].anv[2] = smp[i].drv[2]; + } + + /* Setup the white & black point blend factor, that makes sure the white and black */ + /* points are not displaced. */ + for (i = 0; i < nmpts; i++) { + smp[i].naxbf = comp_naxbf(&opts, smp[i]._sv); +//printf("~1 point %d, comp_lvc = %f, naxbf = %f\n",i,comp_lvc(&opts, smp[i]._sv),smp[i].naxbf); + } + + /* Setup the 3D -> 2D tangent conversion ready for guide vector optimization */ + { + double ta[3] = { 50.0, 0.0, 0.0 }; + double tc[3] = { 0.0, 0.0, 0.0 }; + + for (ix = 0; ix < nmpts; ix++) { + /* Coompute a rotation that brings the target point location to 50,0,0 */ + icmVecRotMat(smp[ix].m2d, smp[ix].sv, sc_gam->cent, ta, tc); + + /* And inverse */ + icmVecRotMat(smp[ix].m3d, ta, tc, smp[ix].sv, sc_gam->cent); + } + } + + /* Figure out which neighbors of the source values to use */ + /* for the relative error calculations. */ + /* Locate the neighbor within the radius for this point, */ + /* and weight them with a Gausian filter weight. */ + /* The radius is computed on the normalised surface for this point. */ + VA(("Establishing filter neighbourhoods\n")); + { + double mm[3][4]; /* Tangent alignment rotation */ + double m2[2][2]; /* Additional matrix to alight a with L axis */ + double ta[3] = { 50.0, 0.0, 0.0 }; + double tc[3] = { 0.0, 0.0, 0.0 }; + double avgnd = 0.0; /* Total the average number of neighbours */ + int minnd = 1e6; /* Minimum number of neighbours */ + + for (ix = 0; ix < nmpts; ix++) { + double tt[3], rrdl, rrdh, rrdc, dd; + double msv[3], ndx[4]; /* Midpoint source value, quadrant distance */ + double pr; /* Average point radius */ + +//printf("~1 computing neigbourhood for point %d at %f %f %f\n",ix, smp[ix].sv[0], smp[ix].sv[1], smp[ix].sv[2]); + /* Compute a rotation that brings the target point location to 50,0,0 */ + icmNormalize33(tt, smp[ix].sv, smp[ix].sgam->cent, 1.0); + icmVecRotMat(mm, tt, smp[ix].sgam->cent, ta, tc); + + /* Add another rotation to orient it so that [1] corresponds */ + /* with the L direction, and [2] corresponds with the */ + /* hue direction. */ + m2[0][0] = m2[1][1] = 1.0; + m2[0][1] = m2[1][0] = 0.0; + tt[0] = smp[ix].sv[0] + 1.0; + tt[1] = smp[ix].sv[1]; + tt[2] = smp[ix].sv[2]; + icmNormalize33(tt, tt, smp[ix].sgam->cent, 1.0); + icmMul3By3x4(tt, mm, tt); + dd = tt[1] * tt[1] + tt[2] * tt[2]; + if (dd > 1e-6) { /* There is a sense of L direction */ + + /* Create the 2x2 rotation matrix */ + dd = sqrt(dd); + tt[1] /= dd; + tt[2] /= dd; + + m2[0][0] = m2[1][1] = tt[1]; + m2[0][1] = tt[2]; + m2[1][0] = -tt[2]; + } + + /* Make rr inversely proportional to radius, so that */ + /* filter scope is constant delta E */ + rrdl = smp[ix].wt.r.rdl; + rrdh = smp[ix].wt.r.rdh; +//printf("~1 rdl %f, rdh %f\n",smp[ix].wt.r.rdl,smp[ix].wt.r.rdh); + if (rrdl < 1e-3) rrdl = 1e-3; + if (rrdh < 1e-3) rrdh = 1e-3; + + /* Average radius of source and destination */ + pr = 0.5 * (smp[ix]._sr + smp[ix].dr); + +//printf("~1 pr = %f from _sr %f & dr %f\n",pr,smp[ix]._sr,smp[ix].dr); + if (pr < 5.0) + pr = 5.0; + rrdl *= 50.0/pr; + rrdh *= 50.0/pr; + rrdc = 0.5 * (rrdl + rrdh); /* Chrominance radius */ + + /* Scale the filter radius by the L/C value, */ + /* so that the filters are largest at the equator, and smallest */ + /* at the white & black points. This allows the wt.a.lx wt.a.cx to work. */ + pr = smp[ix].naxbf + 0.1; /* "spherical" type weighting, 0.707 at 45 degrees */ + rrdl *= pr; + rrdh *= pr; + rrdc *= pr; +//printf("~1 at %f %f %f, rrdl = %f, rrdh = %f\n",smp[ix]._sv[0], smp[ix]._sv[1], smp[ix]._sv[2], rrdl, rrdh); + + smp[ix].nnd = 0; /* Nothing in lists */ + + /* Search for points within the gausian radius */ + for (i = 0; i < nmpts; i++) { + double x, y, z, tv[3]; + + /* compute rotated location */ + icmNormalize33(tt, smp[i].sv, smp[ix].sgam->cent, 1.0); + icmMul3By3x4(tv, mm, tt); + icmMulBy2x2(&tv[1], m2, &tv[1]); + + x = tv[1]/rrdl; + y = tv[2]/rrdh; + z = (tv[0] - 50.0)/rrdc; + + /* Compute normalized delta normalized tangent surface */ + dd = x * x + y * y + z * z; + + /* If we're within the direction filtering radius, */ + /* and not of the opposite hue */ + if (dd <= 1.0 && tv[0] > 0.0) { + double w; + + dd = sqrt(dd); /* Convert to radius <= 1.0 */ + + /* Add this point into the list */ + if (smp[ix].nnd >= smp[ix]._nnd) { + neighb *nd; + int _nnd; + _nnd = 5 + smp[ix]._nnd * 2; + if ((nd = (neighb *)realloc(smp[ix].nd, _nnd * sizeof(neighb))) == NULL) { + VB(("realloc of neighbs at vector %d failed\n",ix)); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + free_nearsmth(smp, nmpts); + *npp = 0; + return NULL; + } + smp[ix].nd = nd; + smp[ix]._nnd = _nnd; + } + smp[ix].nd[smp[ix].nnd].n = &smp[i]; + + /* Box filter */ +// w = 1.0; + + /* Triangle filter */ +// w = 1.0 - dd; + +// /* Cubic spline filter */ +// w = 1.0 - dd; +// w = w * w * (3.0 - 2.0 * w); + + /* Gaussian filter (default) */ + w = exp(-9.0 * dd/2.0); + + /* Sphere filter */ +// w = sqrt(1.0 - dd * dd); + + /* Sinc^2 filter */ +// w = 3.1415926 * dd; +// if (w < 1e-9) +// w = 1e-9; +// w = sin(w)/w; +// w = w * w; + + smp[ix].nd[smp[ix].nnd].w = w; /* Will be normalized to sum to 1.0 */ + +// /* Sphere filter for depth */ +// w = sqrt(1.0 - dd * dd); + + /* Cubic spline filter for depth */ +// w = 1.0 - dd; +// w = w * w * (3.0 - 2.0 * w); + + /* Gaussian filter for depth (default) */ + w = exp(-9.0 * dd/2.0); + + smp[ix].nd[smp[ix].nnd].rw = w; /* Won't be normalized */ + +//printf("~1 adding %d at %f %f %f, rad %f L %f, w %f dir.\n",i, smp[i].sv[0], smp[i].sv[1], smp[i].sv[2],sqrt(dd),tv[0],smp[ix].nd[smp[ix].nnd].w); + + smp[ix].nnd++; + } + } + if (smp[ix].nnd < minnd) + minnd = smp[ix].nnd; + avgnd += (double)smp[ix].nnd; +//printf("~1 total of %d dir neigbours\n\n",smp[ix].nnd); + + } + avgnd /= (double)nmpts; + + if (verb) printf("Average number of direction guide neigbours = %f, min = %d\n",avgnd,minnd); + + /* Now normalize each points weighting */ + for (i = 0; i < nmpts; i++) { + double tw; + + /* Adjust direction weights to sum to 1.0 */ + for (tw = 0.0, j = 0; j < smp[i].nnd; j++) { + tw += smp[i].nd[j].w; + } + for (j = 0; j < smp[i].nnd; j++) { + smp[i].nd[j].w /= tw; + } + + } + } + +#ifdef SHOW_NEIGB_WEIGHTS + { + vrml *wrl = NULL; + double yellow[3] = { 1.0, 1.0, 0.0 }; + double red[3] = { 1.0, 0.0, 0.0 }; + double green[3] = { 0.0, 1.0, 0.0 }; + double pp[3]; + + for (i = 0; i < nmpts; i++) { + double maxw; + + if ((wrl = new_vrml("weights.wrl", 1)) == NULL) + error("New vrml failed"); + + maxw = 0.0; + for (j = 0; j < smp[i].nnd; j++) { + if (smp[i].nd[j].w > maxw) + maxw = smp[i].nd[j].w; + } + for (j = 0; j < smp[i].nnd; j++) { + wrl->add_col_vertex(wrl, 0, smp[i].sgam->cent, smp[i].nd[j].n == &smp[i] ? red : yellow); + icmNormalize33(pp, smp[i].nd[j].n->_sv, smp[i].sgam->cent, smp[i].nd[j].w * 50.0/maxw); + wrl->add_col_vertex(wrl, 0, pp, smp[i].nd[j].n == &smp[i] ? red : yellow); + } + wrl->make_lines(wrl, 0, 2); + + wrl->del(wrl); + printf("Waiting for input after writing 'weights.wrl' for point %d:\n",i); + getchar(); + } + } +#endif /* SHOW_NEIGB_WEIGHTS */ + + /* Optimise the location of the source to destination mapping. */ + if (verb) printf("Optimizing source to destination mapping...\n"); + + VA(("Doing first pass to locate the nearest point\n")); + /* First pass to locate the weighted nearest point, to use in subsequent passes */ + { + double s[2] = { 20.0, 20.0 }; /* 2D search area */ + double iv[3]; /* Initial start value */ + double nv[2]; /* 2D New value */ + double tp[3]; /* Resultint value */ + double ne; /* New error */ + int notrials = NO_TRIALS; + + for (i = 0; i < nmpts; i++) { /* Move all the points */ + double bnv[2]; /* Best 2d value */ + double brv; /* Best return value */ + int trial; + double mv; + + opts.pass = 0; /* Itteration pass */ + opts.ix = i; /* Point to optimise */ + opts.p = &smp[i]; + + /* If the img point is within the destination, then we're */ + /* expanding, so temporarily swap src and radial dest. */ + /* (??? should we use the cvect() direction to determine swap, */ + /* rather than radial ???) */ + smp[i].swap = 0; + if (useexp && smp[i].dr > (smp[i].sr + 1e-9)) { + gamut *tt; + double dd; + + smp[i].swap = 1; + tt = smp[i].dgam; smp[i].dgam = smp[i].sgam; smp[i].sgam = tt; + + smp[i].dr = smp[i].sr; + smp[i].dv[0] = smp[i].sv[0]; + smp[i].dv[1] = smp[i].sv[1]; + smp[i].dv[2] = smp[i].sv[2]; + + smp[i].sr = smp[i].drr; + smp[i].sv[0] = smp[i].drv[0]; + smp[i].sv[1] = smp[i].drv[1]; + smp[i].sv[2] = smp[i].drv[2]; + } + + /* Convert our start value from 3D to 2D for speed. */ + icmMul3By3x4(iv, smp[i].m2d, smp[i].dv); + nv[0] = iv[0] = iv[1]; + nv[1] = iv[1] = iv[2]; + + /* Do several trials from different starting points to avoid */ + /* any local minima, particularly with nearest mapping. */ + brv = 1e38; + for (trial = 0; trial < notrials; trial++) { + double rv; /* Temporary */ + + /* Optimise the point */ + if (powell(&rv, 2, nv, s, 0.01, 1000, optfunc1, (void *)(&opts), NULL, NULL) == 0 + && rv < brv) { + brv = rv; +//printf("~1 point %d, trial %d, new best %f\n",i,trial,sqrt(rv)); + bnv[0] = nv[0]; + bnv[1] = nv[1]; + } +//else printf("~1 powell failed with rv = %f\n",rv); + /* Adjust the starting point with a random offset to avoid local minima */ + nv[0] = iv[0] + d_rand(-20.0, 20.0); + nv[1] = iv[1] + d_rand(-20.0, 20.0); + } + if (brv == 1e38) { /* We failed to get a result */ + VB(("multiple powells failed to get a result\n")); +#ifdef NEVER + /* Optimise the point with debug on */ + opts.debug = 1; + icmMul3By3x4(iv, smp[i].m2d, smp[i].dv); + nv[0] = iv[0] = iv[1]; + nv[1] = iv[1] = iv[2]; + powell(NULL, 2, nv, s, 0.01, 1000, optfunc1, (void *)(&opts), NULL, NULL); +#endif + free_nearsmth(smp, nmpts); + *npp = 0; + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + return NULL; + } + + /* Convert best result 2D -> 3D */ + tp[2] = bnv[1]; + tp[1] = bnv[0]; + tp[0] = 50.0; + icmMul3By3x4(tp, smp[i].m3d, tp); + + /* Remap it to the destinaton gamut surface */ + smp[i].dgam->radial(smp[i].dgam, tp, tp); + icmCpy3(smp[i].aodv, tp); + + /* Undo any swap */ + if (smp[i].swap) { + gamut *tt; + double dd; + + tt = smp[i].dgam; smp[i].dgam = smp[i].sgam; smp[i].sgam = tt; + + /* We get the point on the real src gamut out when swap */ + smp[i]._sv[0] = smp[i].aodv[0]; + smp[i]._sv[1] = smp[i].aodv[1]; + smp[i]._sv[2] = smp[i].aodv[2]; + + /* So we need to compute cusp mapped sv */ + comp_ce(&opts, smp[i].sv, smp[i]._sv, &smp[i].wt); + smp[i].sr = icmNorm33(smp[i].sv, smp[i].sgam->cent); + + VB(("Exp Src %d = %f %f %f\n",i,smp[i]._sv[0],smp[i]._sv[1],smp[i]._sv[2])); + smp[i].aodv[0] = smp[i].drv[0]; + smp[i].aodv[1] = smp[i].drv[1]; + smp[i].aodv[2] = smp[i].drv[2]; + } + } + if (verb) { + printf("."); fflush(stdout); + } + } + + VA(("Locating weighted mapping vectors without smoothing:\n")); + /* Second pass to locate the optimized overall weighted point nrdv[], */ + /* not counting relative error. */ + { + double s[2] = { 20.0, 20.0 }; /* 2D search area */ + double iv[3]; /* Initial start value */ + double nv[2]; /* 2D New value */ + double tp[3]; /* Resultint value */ + double ne; /* New error */ + int notrials = NO_TRIALS; + + for (i = 0; i < nmpts; i++) { /* Move all the points */ + double bnv[2]; /* Best 2d value */ + double brv; /* Best return value */ + int trial; + double mv; + + opts.pass = 0; /* Itteration pass */ + opts.ix = i; /* Point to optimise */ + opts.p = &smp[i]; + +//printf("~1 point %d, sv %f %f %f\n",i,smp[i].sv[0],smp[i].sv[1],smp[i].sv[2]); + + /* Convert our start value from 3D to 2D for speed. */ + icmMul3By3x4(iv, smp[i].m2d, smp[i].aodv); + + nv[0] = iv[0] = iv[1]; + nv[1] = iv[1] = iv[2]; +//printf("~1 point %d, iv %f %f %f, 2D %f %f\n",i, smp[i].aodv[0], smp[i].aodv[1], smp[i].aodv[2], iv[0], iv[1]); + + /* Do several trials from different starting points to avoid */ + /* any local minima, particularly with nearest mapping. */ + brv = 1e38; + for (trial = 0; trial < notrials; trial++) { + double rv; /* Temporary */ + + /* Optimise the point */ + if (powell(&rv, 2, nv, s, 0.01, 1000, optfunc2, (void *)(&opts), NULL, NULL) == 0 + && rv < brv) { + brv = rv; +//printf("~1 point %d, trial %d, new best %f at xy %f %f\n",i,trial,sqrt(rv), nv[0],nv[1]); + bnv[0] = nv[0]; + bnv[1] = nv[1]; + } +//else printf("~1 powell failed with rv = %f\n",rv); + /* Adjust the starting point with a random offset to avoid local minima */ + nv[0] = iv[0] + d_rand(-20.0, 20.0); + nv[1] = iv[1] + d_rand(-20.0, 20.0); + } + if (brv == 1e38) { /* We failed to get a result */ + VB(("multiple powells failed to get a result\n")); +#ifdef NEVER + /* Optimise the point with debug on */ + opts.debug = 1; + icmMul3By3x4(iv, smp[i].m2d, smp[i].dv); + nv[0] = iv[0] = iv[1]; + nv[1] = iv[1] = iv[2]; + powell(NULL, 2, nv, s, 0.01, 1000, optfunc2, (void *)(&opts), NULL, NULL); +#endif + free_nearsmth(smp, nmpts); + *npp = 0; + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + return NULL; + } + + /* Convert best result 3D -> 2D */ + tp[2] = bnv[1]; + tp[1] = bnv[0]; + tp[0] = 50.0; + icmMul3By3x4(tp, smp[i].m3d, tp); + + /* Remap it to the destinaton gamut surface */ + smp[i].dgam->radial(smp[i].dgam, tp, tp); + + icmCpy3(smp[i].nrdv, tp); /* Non smoothed result */ + icmCpy3(smp[i].anv, tp); /* Starting point for smoothing */ + icmCpy3(smp[i].dv, tp); /* Default current solution */ + smp[i].dr = icmNorm33(smp[i].dv, smp[i].dgam->cent); +//printf("~1 %d: dv %f %f %f\n", i, smp[i].dv[0], smp[i].dv[1], smp[i].dv[2]); + + } + if (verb) { + printf("."); fflush(stdout); + } + } + +#ifdef DIAG_POINTS + /* Show just the closest vectors etc. */ + for (i = 0; i < nmpts; i++) { /* Move all the points */ +// icmCpy3(smp[i].dv, smp[i].drv); /* Radial */ + icmCpy3(smp[i].dv, smp[i].aodv); /* Nearest */ +// icmCpy3(smp[i].dv, smp[i].nrdv); /* No smoothed weighted */ +// icmCpy3(smp[i].dv, smp[i].dv); /* pre-filter smooothed */ + smp[i].dr = icmNorm33(smp[i].dv, smp[i].dgam->cent); + } +#else + /* The smoothed direction and raw depth is a single pass, */ + /* but we use multiple passes to determine the extra depth that */ + /* needs to be added so that the smoothed result lies within */ + /* the destination gamut. */ + +#if VECADJPASSES > 0 || RSPLPASSES > 0 + /* We will need inward pointing correction vectors */ + { + gamut *shgam; /* Shrunken di_gam */ + double cusps[6][3]; + double wp[3], bp[3], kp[3]; + double p[3], p2[3], rad; + int i; + + double s[2] = { 20.0, 20.0 }; /* 2D search area */ + double iv[3]; /* Initial start value */ + double nv[2]; /* 2D New value */ + double tp[3]; /* Resultint value */ + double ne; /* New error */ + int notrials = NO_TRIALS; + + cow *gpnts = NULL; /* Mapping points to create 3D -> 3D mapping */ + datai il, ih; + datao ol, oh; + int gres[MXDI]; + double avgdev[MXDO]; + + /* Create a gamut that is a shrunk version of the destination */ + + if ((shgam = new_gamut(di_gam->getsres(di_gam), di_gam->getisjab(di_gam), + di_gam->getisrast(di_gam))) == NULL) { + free_nearsmth(smp, nmpts); + *npp = 0; + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + return NULL; + } + + shgam->setnofilt(shgam); + + /* Translate all the surface nodes */ + for (i = 0;;) { + double len; + + if ((i = di_gam->getrawvert(di_gam, p, i)) < 0) + break; + + doshrink(&opts, p, p, SHRINK); + shgam->expand(shgam, p); + } + /* Translate cusps */ + if (di_gam->getcusps(di_gam, cusps) == 0) { + shgam->setcusps(shgam, 0, NULL); + for (i = 0; i < 6; i++) { + doshrink(&opts, p, cusps[i], SHRINK); + shgam->setcusps(shgam, 1, p); + } + shgam->setcusps(shgam, 2, NULL); + } + /* Translate white and black points */ + if (di_gam->getwb(di_gam, wp, bp, kp, NULL, NULL, NULL) == 0) { + doshrink(&opts, wp, wp, SHRINK); + doshrink(&opts, bp, bp, SHRINK); + doshrink(&opts, kp, kp, SHRINK); + shgam->setwb(shgam, wp, bp, kp); + } + + if ((gpnts = (cow *)malloc(nmpts * sizeof(cow))) == NULL) { + fprintf(stderr,"gamut map: Malloc of near smooth points failed\n"); + free_nearsmth(smp, nmpts); + *npp = 0; + shgam->del(shgam); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + return NULL; + } + + /* Now locate the closest points on the shrunken gamut */ + /* and set them up for creating a rspl */ + opts.shgam = shgam; + for (i = 0; i < nmpts; i++) { /* Move all the points */ + gtri *ctri = NULL; + double tmp[3]; + double bnv[2]; /* Best 2d value */ + double brv; /* Best return value */ + int trial; + double mv; + + opts.pass = 0; /* Itteration pass */ + opts.ix = i; /* Point to optimise */ + opts.p = &smp[i]; + + /* Convert our start value from 3D to 2D for speed. */ + icmMul3By3x4(iv, smp[i].m2d, smp[i].dv); + nv[0] = iv[0] = iv[1]; + nv[1] = iv[1] = iv[2]; + + /* Do several trials from different starting points to avoid */ + /* any local minima, particularly with nearest mapping. */ + brv = 1e38; + for (trial = 0; trial < notrials; trial++) { + double rv; /* Temporary */ + + /* Optimise the point */ + if (powell(&rv, 2, nv, s, 0.01, 1000, optfunc1a, (void *)(&opts), NULL, NULL) == 0 + && rv < brv) { + brv = rv; + bnv[0] = nv[0]; + bnv[1] = nv[1]; + } + /* Adjust the starting point with a random offset to avoid local minima */ + nv[0] = iv[0] + d_rand(-20.0, 20.0); + nv[1] = iv[1] + d_rand(-20.0, 20.0); + } + if (brv == 1e38) { /* We failed to get a result */ + VB(("multiple powells failed to get a result\n")); + shgam->del(shgam); /* Done with this */ + free_nearsmth(smp, nmpts); + *npp = 0; + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + return NULL; + } + + /* Convert best result 2D -> 3D */ + tp[2] = bnv[1]; + tp[1] = bnv[0]; + tp[0] = 50.0; + icmMul3By3x4(tp, smp[i].m3d, tp); + + /* Remap it to the destinaton gamut surface */ + shgam->radial(shgam, tp, tp); + + /* Compute mapping vector from dst to shdst */ + icmSub3(smp[i].temp, tp, smp[i].dv); + + + /* In case shrunk vector is very short, add a small part */ + /* of the nearest normal. */ + smp[i].dgam->nearest_tri(smp[i].dgam, NULL, smp[i].dv, &ctri); + icmScale3(tmp, ctri->pe, 0.1); /* Scale to small inwards */ + icmAdd3(smp[i].temp, smp[i].temp, tmp); + + /* evector */ + icmNormalize3(smp[i].temp, smp[i].temp, 1.0); + + /* Place it in rspl setup array */ + icmCpy3(gpnts[i].p, smp[i].dv); + icmCpy3(gpnts[i].v, smp[i].temp); + gpnts[i].w = 1.0; + } + + for (j = 0; j < 3; j++) { /* Set resolution for all axes */ +// gres[j] = (mapres+1)/2; /* Half resolution */ + gres[j] = mapres; /* Full resolution */ + avgdev[j] = GAMMAP_RSPLAVGDEV; + } + + evectmap = new_rspl(RSPL_NOFLAGS, 3, 3); /* Allocate 3D -> 3D */ + evectmap->fit_rspl_w(evectmap, GAMMAP_RSPLFLAGS, gpnts, nmpts, + map_il, map_ih, gres, map_ol, map_oh, 1.0, avgdev, NULL); + +// ~~999 +#ifdef PLOT_EVECTS /* Create VRML of error correction vectors */ + { + vrml *wrl = NULL; + int doaxes = 0; + double cc[3] = { 0.7, 0.7, 0.7 }; + double red[3] = { 1.0, 0.0, 0.0 }; + double green[3] = { 0.0, 1.0, 0.0 }; + double tmp[3]; + co cp; +#ifdef PLOT_AXES + doaxes = 1; +#endif + + printf("###### gamut/nearsmth.c: writing diagnostic evects.wrl\n"); + wrl = new_vrml("evects.wrl", doaxes); + wrl->make_gamut_surface_2(wrl, di_gam, 0.6, 0, cc); + cc[0] = -1.0; + wrl->make_gamut_surface(wrl, shgam, 0.2, cc); + + /* Start of guide vector plot */ + wrl->start_line_set(wrl, 0); + + for (i = 0; i < nmpts; i++) { + wrl->add_col_vertex(wrl, 0, smp[i].dv, red); +#ifdef NEVER /* Plot created vectors */ + icmScale3(tmp, smp[i].temp, 4.0); + icmSub3(tmp, smp[i].dv, tmp); +#else + /* Plot interpolated vectors */ + icmCpy3(cp.p, smp[i].dv); + evectmap->interp(evectmap, &cp); + icmScale3(tmp, cp.v, 4.0); + icmSub3(tmp, smp[i].dv, tmp); +#endif + wrl->add_col_vertex(wrl, 0, tmp, green); + } + wrl->make_lines(wrl, 0, 2); /* Guide vectors */ + wrl->del(wrl); /* Write and delete */ + } +#endif /* PLOT_EVECTS */ + shgam->del(shgam); /* Done with this */ + free(gpnts); + } +#endif /* VECADJPASSES > 0 || RSPLPASSES > 0 */ + +#ifdef VECSMOOTHING + VA(("Smoothing guide vectors:\n")); + + /* Compute the neighbourhood smoothed anv[] from dv[] */ + for (i = 0; i < nmpts; i++) { + double anv[3]; /* new anv[] */ + double tmp[3]; + + /* Compute filtered value */ + anv[0] = anv[1] = anv[2] = 0.0; + for (j = 0; j < smp[i].nnd; j++) { + nearsmth *np = smp[i].nd[j].n; /* Pointer to neighbor */ + double nw = smp[i].nd[j].w; /* Weight */ + double tmp[3]; + + icmSub3(tmp, smp[i].sv, np->sv); /* Vector from neighbour src to src */ + icmAdd3(tmp, tmp, np->dv); /* Neigbour dst + vector */ + + icmScale3(tmp, tmp, nw); /* weight for filter */ + icmAdd3(anv, anv, tmp); /* sum filtered value */ + } + + /* Blend to un-smoothed value on neutral axis */ + icmBlend3(anv, smp[i].dv, anv, smp[i].naxbf); + + icmCpy3(smp[i].dv, anv); + icmCpy3(smp[i].anv, anv); + smp[i].rext = 0.0; /* No correction */ + } + +#if VECADJPASSES > 0 + /* Fine tune vectors to compensate for side effects of vector smoothing */ + + VA(("Fine tuning out of gamut guide vectors:\n")); + + /* Loopkup correction vectors */ + VA(("Computing fine tuning direction:\n")); + for (i = 0; i < nmpts; i++) { + co cp; + double nd, id, tmp[3]; + + icmCpy3(cp.p, smp[i].dv); + evectmap->interp(evectmap, &cp); + icmNormalize3(smp[i].evect, cp.v, 1.0); + + /* ~~99 ?? should we deal with white & black direction here ?? */ + + /* Use closest as a default */ + smp[i].dgam->nearest(smp[i].dgam, smp[i].tdst, smp[i].dv); + nd = icmNorm33(smp[i].tdst, smp[i].dv); /* Dist to nearest */ + + /* Compute intersection with dest gamut as tdst */ + if (!vintersect2(smp[i].dgam, NULL, tmp, smp[i].evect, smp[i].dv)) { + /* Got an intersection */ + id = icmNorm33(tmp, smp[i].dv); /* Dist to intersection */ + if (id <= (nd + 5.0)) /* And it seems sane */ + icmCpy3(smp[i].tdst, tmp); + } + + smp[i].rext = 0.0; + } + + VA(("Fine tuning guide vectors:\n")); + for (it = 0; it < VECADJPASSES; it++) { + double avgog = 0.0, maxog = 0.0, nog = 0.0; + double avgig = 0.0, maxig = 0.0, nig = 0.0; + + /* Filter the level of out/in gamut, and apply correction vector */ + for (i = 0; i < nmpts; i++) { + double cvec[3], clen; + double minext = 1e80; + double maxext = -1e80; /* Max weighted depth extension */ + double dext, gain; + + minext = -20.0; + + /* Compute filtered value */ + for (j = 0; j < smp[i].nnd; j++) { + nearsmth *np = smp[i].nd[j].n; /* Pointer to neighbor */ + double nw = smp[i].nd[j].rw; /* Weight */ + double tmpl; + + icmSub3(cvec, np->tdst, np->anv); /* Vector needed to target for neighbour */ + clen = icmDot3(smp[i].evect, cvec); /* Error in this direction */ + + tmpl = nw * (clen - minext); /* Track maximum weighted extra depth */ + if (tmpl < 0.0) + tmpl = 0.0; + if (tmpl > maxext) + maxext = tmpl; + } + maxext += minext; + + if (it == 0) + gain = 1.2; + else + gain = 0.8; + + /* Accumulate correction with damping */ + smp[i].rext += gain * maxext; + + /* Error for just this point */ + icmSub3(cvec, smp[i].tdst, smp[i].anv); + clen = icmDot3(smp[i].evect, cvec); + + /* Blend to individual correction on neutral axis */ + dext = smp[i].naxbf * smp[i].rext + (1.0 - smp[i].naxbf) * clen; + + /* Apply integrated correction */ + icmScale3(cvec, smp[i].evect, dext); + icmAdd3(smp[i].anv, smp[i].dv, cvec); + + if (clen > 0.0) { /* Compression */ + if (clen > maxog) + maxog = clen; + avgog += clen; + nog++; + + } else { /* Expansion */ + if (-clen > maxig) + maxig = -clen; + avgig += -clen; + nig++; + } + } + if (verb) + printf("No og %4.0f max %f avg %f, No ig %4.0f max %f avg %f\n", + nog,maxog,nog > 1 ? avgog/nog : 0.0, nig,maxig,nig > 1 ? avgig/nig : 0.0); + } + + /* Copy final results */ + for (i = 0; i < nmpts; i++) { + icmCpy3(smp[i].dv, smp[i].anv); + smp[i].dr = icmNorm33(smp[i].dv, smp[i].dgam->cent); + } + + if (verb) { + double avgog = 0.0, maxog = 0.0, nog = 0.0; + double avgig = 0.0, maxig = 0.0, nig = 0.0; + + /* Check the result */ + for (i = 0; i < nmpts; i++) { + double cvec[3], clen; + + /* Error for just this point, for stats */ + icmSub3(cvec, smp[i].tdst, smp[i].anv); + clen = icmDot3(smp[i].evect, cvec); + + if (clen > 0.0) { /* Compression */ + if (clen > maxog) + maxog = clen; + avgog += clen; + nog++; + + } else { /* Expansion */ + if (-clen > maxig) + maxig = -clen; + avgig += -clen; + nig++; + } + } + printf("No og %4.0f max %f avg %f, No ig %4.0f max %f avg %f\n", + nog,maxog,nog > 1 ? avgog/nog : 0.0, nig,maxig,nig > 1 ? avgig/nig : 0.0); + } +#endif /* VECADJUST */ +#endif /* VECADJPASSES > 0 */ + +#if RSPLPASSES > 0 + /* We need to adjust the vectors with extra depth to compensate for */ + /* for the effect of rspl smoothing. */ + { + cow *gpnts = NULL; /* Mapping points to create 3D -> 3D mapping */ + rspl *map = NULL; /* Test map */ + datai il, ih; + datao ol, oh; + int gres[MXDI]; + double avgdev[MXDO]; + double icgain, ixgain; /* Initial compression, expansion gain */ + double fcgain, fxgain; /* Final compression, expansion gain */ + + VA(("Fine tuning vectors to allow for rspl smoothing:\n")); + + for (j = 0; j < 3; j++) { /* Copy ranges */ + il[j] = map_il[j]; + ih[j] = map_ih[j]; + ol[j] = map_ol[j]; + oh[j] = map_oh[j]; + } + + /* Adjust the input ranges for guide vectors */ + for (i = 0; i < nmpts; i++) { + for (j = 0; j < 3; j++) { + if (smp[i]._sv[j] < il[j]) + il[j] = smp[i]._sv[j]; + if (smp[i]._sv[j] > ih[j]) + ih[j] = smp[i]._sv[j]; + } + } + + /* Now expand the bounding box by aprox 5% margin, but scale grid res */ + /* to match, so that the natural or given boundary still lies on the grid. */ + /* (This duplicates code in gammap applied after near_smooth() returns) */ + /* (We are assuming that our changes to the giude vectprs won't expand the ranges) */ + { + int xmapres; + double scale; + + xmapres = (int) ((mapres-1) * 0.05 + 0.5); + if (xmapres < 1) + xmapres = 1; + + scale = (double)(mapres-1 + xmapres)/(double)(mapres-1); + + for (j = 0; j < 3; j++) { + double low, high; + high = ih[j]; + low = il[j]; + ih[j] = (scale * (high - low)) + low; + il[j] = (scale * (low - high)) + high; + } + + mapres += 2 * xmapres; + } + + if ((gpnts = (cow *)malloc(nmpts * sizeof(cow))) == NULL) { + fprintf(stderr,"gamut map: Malloc of near smooth points failed\n"); + free_nearsmth(smp, nmpts); + *npp = 0; + if (evectmap != NULL) + evectmap->del(evectmap); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + return NULL; + } + + /* Loopkup correction vectors */ + VA(("Computing fine tuning direction for vectors:\n")); + for (i = 0; i < nmpts; i++) { + co cp; + double nd, id, tmp[3]; + + icmCpy3(cp.p, smp[i].dv); + evectmap->interp(evectmap, &cp); + icmNormalize3(smp[i].evect, cp.v, 1.0); + + /* ~~99 ?? should we deal with white & black direction here ?? */ + + /* Use closest as a default */ + smp[i].dgam->nearest(smp[i].dgam, smp[i].tdst, smp[i].dv); + nd = icmNorm33(smp[i].tdst, smp[i].dv); /* Dist to nearest */ + + /* Compute intersection with dest gamut as tdst */ + if (!vintersect2(smp[i].dgam, NULL, tmp, smp[i].evect, smp[i].dv)) { + /* Got an intersection */ + id = icmNorm33(tmp, smp[i].dv); /* Dist to intersection */ + if (id <= (nd + 5.0)) /* And it seems sane */ + icmCpy3(smp[i].tdst, tmp); + } + + smp[i].coff[0] = smp[i].coff[1] = smp[i].coff[2] = 0.0; + smp[i].rext = 0.0; + } + + /* We know initially that dv == anv */ + /* Each pass computes a rext for each point, then */ + /* anv[] = dv[] + smooth(rext * evect[]) to try and avoid clipping */ + VA(("Fine tune guide vectors for rspl:\n")); + for (it = 0; it < RSPLPASSES; it++) { + double tmp[3]; + double avgog = 0.0, maxog = 0.0, nog = 0.0; + double avgig = 0.0, maxig = 0.0, nig = 0.0; + double avgrext = 0.0; + double ovlen; + + /* Setup the rspl guide points for creating rspl */ + for (i = 0; i < nmpts; i++) { + icmCpy3(gpnts[i].p, smp[i]._sv); /* The orgininal src point */ + icmCpy3(gpnts[i].v, smp[i].anv); /* current dst from previous results */ + gpnts[i].w = 1.0; + } + + for (j = 0; j < 3; j++) { /* Set resolution for all axes */ +// gres[j] = (mapres+1)/2; /* Half resolution for speed */ + gres[j] = mapres; /* Full resolution */ + avgdev[j] = GAMMAP_RSPLAVGDEV; + } + map = new_rspl(RSPL_NOFLAGS, 3, 3); /* Allocate 3D -> 3D */ + map->fit_rspl_w(map, GAMMAP_RSPLFLAGS, gpnts, nmpts, + il, ih, gres, ol, oh, mapsmooth, avgdev, NULL); + + /* See what the source actually maps to via rspl, and how far from */ + /* the target point they are. */ + for (i = 0; i < nmpts; i++) { + co cp; + double cvec[3]; + + /* Lookup rspl smoothed destination value */ + icmCpy3(cp.p, smp[i]._sv); + map->interp(map, &cp); + icmCpy3(smp[i].temp, cp.v); + + /* Compute the correction needed and it's signed length */ + icmSub3(cvec, smp[i].tdst, smp[i].temp); + smp[i].clen = icmDot3(smp[i].evect, cvec); + } + + /* Compute local correction */ + for (i = 0; i < nmpts; i++) { + double minext = 1e80; + double maxext = -1e80; /* Max weighted depth extension */ + double clen; + double tpoint[3], cvect[3]; + double tt; + double cgain, xgain; /* This itters compression, expansion gain */ + double gain; /* Gain used */ + + /* See what the worst case is in the local area, and */ + /* aim to lower the whole local area by enough to */ + /* cause the max to be 0.0 (just on the gamut) */ + + /* Compute local depth value */ + minext = -20.0; /* Base to measure max from */ + for (j = 0; j < smp[i].nnd; j++) { + nearsmth *np = smp[i].nd[j].n; /* Pointer to neighbor */ + double nw = smp[i].nd[j].rw; /* Weight */ + double tmpl; + + tmpl = nw * (np->clen - minext); /* Track maximum weighted extra depth */ + if (tmpl < 0.0) + tmpl = 0.0; + if (tmpl > maxext) + maxext = tmpl; + } + maxext += minext; + + /* maxext is the current effective error at this point with rext aim point */ + if (it == 0) { /* Set target on first itteration */ + if (smp[i].rext <= 0.0) /* Expand direction */ + smp[i].rext += maxext; + else + smp[i].rext += RSPLSCALE * maxext; + } + avgrext += smp[i].rext; + + /* Compute offset target point at maxlen from tdst in evect dir. */ + icmScale3(tpoint, smp[i].evect, smp[i].rext); + icmAdd3(tpoint, tpoint, smp[i].tdst); + + /* Expansion/compression gain program */ + icgain = 1.4; /* Initial itteration compression gain */ + ixgain = smp[i].wt.f.x * icgain; /* Initial itteration expansion gain */ + fcgain = 0.5 * icgain; /* Final itteration compression gain */ + fxgain = 0.5 * ixgain; /* Final itteration expansion gain */ + + /* Set the gain */ + tt = it/(RSPLPASSES - 1.0); + cgain = (1.0 - tt) * icgain + tt * fcgain; + xgain = (1.0 - tt) * ixgain + tt * fxgain; + if (it != 0) /* Expand only on first itter */ + xgain = 0.0; + +//if (i == 0) printf("~1 i %d, it %d, wt.f.x = %f, cgain %f, xgain %f\n",i,it,smp[i].wt.f.x, cgain, xgain); + if (smp[i].rext > 0.0) /* Compress direction */ + gain = cgain; + else + gain = xgain; + + /* Keep stats of this point */ + clen = smp[i].clen; + if (clen > 0.0) { + if (clen > maxog) + maxog = clen; + avgog += clen; + nog++; + + } else { /* Expand */ + if (-clen > maxig) + maxig = -clen; + avgig += -clen; + nig++; + } + + /* Compute needed correction from current rspl smoothed anv */ + /* to offset target point. */ + icmSub3(cvect, tpoint, smp[i].temp); /* Correction still needed */ + icmScale3(cvect, cvect, gain); /* Times gain */ + icmAdd3(smp[i].coff, smp[i].coff, cvect); /* Accumulated */ + + icmCpy3(gpnts[i].p, smp[i].dv); + icmCpy3(gpnts[i].v, smp[i].coff); + gpnts[i].w = 1.0; + } + + for (j = 0; j < 3; j++) { /* Set resolution for all axes */ +// gres[j] = (mapres+1)/2; /* Half resolution */ + gres[j] = mapres; /* Full resolution */ + avgdev[j] = GAMMAP_RSPLAVGDEV; + } + map = new_rspl(RSPL_NOFLAGS, 3, 3); /* Allocate 3D -> 3D */ + map->fit_rspl_w(map, GAMMAP_RSPLFLAGS, gpnts, nmpts, + il, ih, gres, ol, oh, 2.0, avgdev, NULL); + + /* Lookup the smoothed extension vector for each point and apply it */ + for (i = 0; i < nmpts; i++) { + double tt; + co cp; + + icmCpy3(cp.p, smp[i].dv); + map->interp(map, &cp); +#ifdef RSPLUSEPOW + spow3(smp[i].coff, cp.v, 1.0/2.0); /* Filtered value is current value */ +#else + icmCpy3(smp[i].coff, cp.v); /* Filtered value is current value */ +#endif + + /* Make sure anv[] is on the destination gamut at neutral axis */ + icmScale3(cp.v, cp.v, smp[i].naxbf); + + /* Apply accumulated offset */ + icmAdd3(smp[i].anv, smp[i].dv, cp.v); + } + map->del(map); + + if (verb) + printf("No og %4.0f max %f avg %f, No ig %4.0f max %f avg %f, avg rext %f\n", + nog,maxog,nog > 1 ? avgog/nog : 0.0, nig,maxig,nig > 1 ? avgig/nig : 0.0, avgrext/nmpts); + } /* Next pass */ + free(gpnts); + + /* Copy last anv to dv for result */ + for (i = 0; i < nmpts; i++) { +// ~~99 +// Normal target + icmCpy3(smp[i].dv, smp[i].anv); + +// Show evect direction +// icmCpy3(smp[i]._sv, smp[i].dv); +// icmScale3(smp[i].evect, smp[i].evect, 4.0); +// icmAdd3(smp[i].dv, smp[i].evect, smp[i].dv); + +// Show target point destination +// icmCpy3(smp[i].dv, smp[i].tdst); + +// Show offset target destination +// icmScale3(smp[i].dv, smp[i].evect, smp[i].rext); +// icmAdd3(smp[i].dv, smp[i].dv, smp[i].tdst); + + smp[i].dr = icmNorm33(smp[i].dv, smp[i].dgam->cent); + } + } + +#endif /* RSPLPASSES > 0 */ +#endif /* NEVER (show debug values) */ + + VA(("Smoothing passes done, doing final houskeeping\n")); + + if (verb) + printf("\n"); + +#if defined(SAVE_VRMLS) && defined(PLOT_MAPPING_INFLUENCE) + create_influence_plot(smp, nmpts); +#endif + + VB(("Final guide points:\n")); + + /* Restore the actual non elevated and cust rotated source point */ + for (i = 0; i < nmpts; i++) { + + VB(("Src %d = %f %f %f\n",i,smp[i].sv[0],smp[i].sv[1],smp[i].sv[2])); + VB(("Dst %d = %f %f %f\n",i,smp[i].dv[0],smp[i].dv[1],smp[i].dv[2])); + + /* Save the cusp mapped source value */ + icmCpy3(smp[i].csv, smp[i].sv); + + /* Finally un cusp map the source point */ +// inv_comp_ce(&opts, smp[i].sv, smp[i].sv, &smp[i].wt); +// smp[i].sr = icmNorm33(smp[i].sv, smp[i].sgam->cent); + icmCpy3(smp[i].sv, smp[i]._sv); + smp[i].sr = smp[i]._sr; + } + + VB(("Creating sub-surface guide points:\n")); + + /* Create sub-surface points. */ + for (i = 0; i < nmpts; i++) { + + /* Create a sub-surface mapping point too. */ + /* Note that not every mapping point has a sub-surface point, */ + /* and that the gflag and vflag will be nz if it does. */ + /* We're assuming here that the dv is close to being on the */ + /* destination gamut, so that the vector_isect param will be */ + /* close to 1.0 at the intended destination gamut. */ + { + double mv[3], ml, nv[3]; /* Mapping vector & length, noralized mv */ + double minv[3], maxv[3]; + double mint, maxt; + gtri *mintri, *maxtri; + + smp[i].vflag = smp[i].gflag = 0; /* Default unknown */ + smp[i].w2 = 0.0; + icmSub3(mv, smp[i].dv, smp[i].sv); /* Mapping vector */ + ml = icmNorm3(mv); /* It's length */ + if (ml > 0.1) { /* If mapping is non trivial */ + +//#define PFCOND i == 802 + +//if (PFCOND) printf("~1 mapping %d = %f %f %f -> %f %f %f\n", i, smp[i].sv[0],smp[i].sv[1],smp[i].sv[2],smp[i].dv[0],smp[i].dv[1],smp[i].dv[2]); +//if (PFCOND) printf("~1 vector %f %f %f, len %f\n", mv[0], mv[1], mv[2],ml); + /* Compute actual depth of ray into destination gamut */ + if (di_gam->vector_isect(di_gam, smp[i].sv, smp[i].dv, + minv, maxv, &mint, &maxt, &mintri, &maxtri) != 0) { + double wp[3], bp[3]; /* Gamut white and black points */ + double p1, napoint[3] = { 50.0, 0.0, 0.0 }; /* Neutral axis point */ + double natarg[3]; /* Neutral axis sub target */ + double adepth1, adepth2 = 1000.0; /* Directional depth, radial depth */ + double adepth; /* Minimum available depth */ + double mv2[3], sml; /* Sub-surface mapping vector & norm. length */ + + /* Locate the point on the neutral axis that is closest to */ + /* the guide ray. We use this as a destination direction */ + /* if the sub surface ray gets very long, and to compute */ + /* a sanity check on the available depth. */ + if (d_gam->getwb(d_gam, NULL, NULL, NULL, wp, dst_kbp ? NULL : bp, dst_kbp ? bp : NULL) == 0) { + if (icmLineLineClosest(napoint, NULL, &p1, NULL, bp, wp, + smp[i].sv,smp[i].dv) == 0) { + /* Clip it */ + if (p1 < 0.0) + icmCpy3(napoint, bp); + else if (p1 > 1.0) + icmCpy3(napoint, wp); + +//if (PFCOND) printf("~1 neutral axis point = %f %f %f\n", napoint[0], napoint[1], napoint[2]); + /* Compute a normalized available depth from distance */ + /* to closest to neautral axis point */ + if (maxt > 1.0) /* Compression */ + if ((mint > 1e-8 && maxt > -1e-8) /* G. & V. Compression */ + || ((mint < -1e-8 && maxt > -1e-8) /* G. Exp & V. comp. */ + && (fabs(mint) < (fabs(maxt) - 1e-8)))) + adepth2 = icmNorm33(napoint, smp[i].dv); + else /* Expansion */ + adepth2 = icmNorm33(napoint, smp[i].sv); + } +#ifdef VERB + else { + printf("icmLineLineClosest failed\n"); + } +#endif + } +#ifdef VERB + else { + printf("d_gam->getwb failed\n"); + } +#endif + +//printf("\n~1 i %d: %f %f %f -> %f %f %f\n isect at t %f and %f\n", i, smp[i].sv[0], smp[i].sv[1], smp[i].sv[2], smp[i].dv[0], smp[i].dv[1], smp[i].dv[2], mint, maxt); + + /* Only create sub-surface mapping vectors if it makes sense. */ + /* If mapping vector is pointing away from destination gamut, */ + /* (which shouldn't happen), ignore it. If the directional depth */ + /* is very thin compared to the radial depth, indicating that we're */ + /* near a "lip", ignore it. */ + if (mint >= -1e-8 && maxt > 1e-8) { + + /* Gamut compression and vector compression */ + if (fabs(mint - 1.0) < fabs(maxt) - 1.0 + && smp[i].dgam->radial(smp[i].dgam, NULL, smp[i].dv) + < smp[i].sgam->radial(smp[i].sgam, NULL, smp[i].dv)) { + +//if (PFCOND) printf("~1 point is gamut comp & vect comp.\n"); +//if (PFCOND) printf("~1 point is gamut comp & vect comp. mint %f maxt %f\n",mint,maxt); + adepth1 = ml * 0.5 * (maxt + mint - 2.0); +#ifdef CYLIN_SUBVEC + adepth = adepth2; /* Always cylindrical depth */ +#else + adepth = adepth1 < adepth2 ? adepth1 : adepth2; /* Smaller of the two */ +#endif + if (adepth1 < (0.5 * adepth2)) + continue; +//if (PFCOND) printf("~1 dir adepth %f, radial adapeth %f\n",adepth1,adepth2); + adepth *= 0.9; /* Can't use 100% */ + smp[i].gflag = 1; /* Gamut compression and */ + smp[i].vflag = 1; /* vector compression */ + + /* Compute available depth and knee factor adjusted sub-vector */ + icmCpy3(smp[i].sv2, smp[i].dv); /* Sub source is guide dest */ + ml *= (1.0 - gamcknf); /* Scale by knee */ + adepth *= (1.0 - gamcknf); + sml = ml < adepth ? ml : adepth; /* Smaller of two */ +//if (PFCOND) printf("~1 adjusted subvec len %f\n",sml); + icmNormalize3(mv2, mv, sml); /* Full sub-surf disp. == no knee */ + icmAdd3(mv2, smp[i].sv2, mv2); /* Knee adjusted destination */ + +//if (PFCOND) printf("~1 before blend sv2 %f %f %f, dv2 %f %f %f\n", smp[i].sv2[0], smp[i].sv2[1], smp[i].sv2[2], mv2[0], mv2[1], mv2[2]); + /* Blend towards n.axis as length of sub vector approaches */ + /* distance to neutral axis. */ + icmSub3(natarg, napoint, smp[i].sv2); + icmNormalize3(natarg, natarg, sml); /* Sub vector towards n.axis */ + icmAdd3(natarg, natarg, smp[i].sv2); /* n.axis target */ +#ifdef CYLIN_SUBVEC + icmCpy3(mv2, natarg); /* cylindrical direction vector */ +#else + icmBlend3(mv2, mv2, natarg, sml/adepth2); +#endif /* CYLIN_SUBVEC */ +//if (PFCOND) printf("~1 after blend sv2 %f %f %f, dv2 %f %f %f\n", smp[i].sv2[0], smp[i].sv2[1], smp[i].sv2[2], mv2[0], mv2[1], mv2[2]); + + icmCpy3(smp[i].dv2, mv2); /* Destination */ + icmCpy3(smp[i].temp, smp[i].dv2); /* Save a copy to temp */ + smp[i].w2 = 0.8; + } else { +//if (PFCOND) printf("~1 point is gamut exp & vect exp. mint %f maxt %f\n",mint,maxt); + smp[i].gflag = 2; /* Gamut expansion and */ + smp[i].vflag = 0; /* vector expansion, */ + /* but crossing over, so no sub vect. */ +//if (PFCOND) printf("~1 point is crossover point\n",mint,maxt); + } + + /* Gamut expansion and vector expansion */ + } else if (mint < -1e-8 && maxt > 1e-8) { + +//if (PFCOND) printf("~1 point is gamut exp & vect exp. mint %f maxt %f\n",mint,maxt); + /* This expand/expand case has reversed src/dst sense to above */ + adepth1 = ml * 0.5 * -mint; +#ifdef CYLIN_SUBVEC + adepth = adepth2; /* Always cylindrical depth */ +#else + adepth = adepth1 < adepth2 ? adepth1 : adepth2; +#endif +//if (PFCOND) printf("~1 dir adepth %f, radial adapeth %f\n",adepth1,adepth2); + adepth *= 0.9; /* Can't use 100% */ + + if (adepth1 < (0.6 * adepth2)) + continue; + + smp[i].gflag = 2; /* Gamut expansion */ + smp[i].vflag = 2; /* vector is expanding */ + + icmCpy3(smp[i].dv2, smp[i].sv); /* Sub dest is guide src */ + ml *= (1.0 - gamxknf); /* Scale by knee */ + adepth *= (1.0 - gamxknf); + sml = ml < adepth ? ml : adepth;/* Smaller of two */ + icmNormalize3(mv2, mv, sml); /* Full sub-surf disp. == no knee */ + icmSub3(mv2, smp[i].dv2, mv2); /* Knee adjusted source */ + + /* Blend towards n.axis as length of sub vector approaches */ + /* distance to neutral axis. */ + icmSub3(natarg, smp[i].dv2, napoint); + icmNormalize3(natarg, natarg, sml); /* Sub vector away n.axis */ + icmSub3(natarg, smp[i].dv2, natarg);/* n.axis oriented source */ +#ifdef CYLIN_SUBVEC + icmCpy3(mv2, natarg); /* cylindrical direction vector */ +#else + icmBlend3(mv2, mv2, natarg, sml/adepth2); /* dir adjusted src */ +#endif /* CYLIN_SUBVEC */ + + icmCpy3(smp[i].sv2, mv2); /* Source */ + icmCpy3(smp[i].temp, smp[i].dv2); /* Save a copy to temp */ + smp[i].w2 = 0.8; + + /* Conflicted case */ + } else { + /* Nonsense vector */ + smp[i].gflag = 0; /* Gamut compression but */ + smp[i].vflag = 0; /* vector is expanding */ +//if (PFCOND) printf("~1 point is nonsense vector mint %f maxt %f\n",mint,maxt); + + icmCpy3(smp[i].dv, smp[i].aodv); /* Clip to the destination gamut */ + } + } + } + } + +#ifdef NEVER // Diagnostic + smp[i].vflag = 0; /* Disable sub-points */ +#endif /* NEVER */ + + VB(("Out Src %d = %f %f %f\n",i,smp[i].sv[0],smp[i].sv[1],smp[i].sv[2])); + VB(("Out Dst %d = %f %f %f\n",i,smp[i].dv[0],smp[i].dv[1],smp[i].dv[2])); + if (smp[i].vflag != 0) { + VB(("Out Src2 %d = %f %f %f\n",i,smp[i].sv2[0],smp[i].sv2[1],smp[i].sv2[2])); + VB(("Out Dst2 %d = %f %f %f\n",i,smp[i].dv2[0],smp[i].dv2[1],smp[i].dv2[2])); + } + } + +#ifdef SUBVEC_SMOOTHING + VB(("Smoothing sub-surface guide points:\n")); + + /* Smooth the sub-surface mapping points */ + /* dv2[] is duplicated in temp[], so use temp[] as the values to be filtered */ + for (i = 0; i < nmpts; i++) { + double tmp[3]; + double fdv2[3]; /* Filtered dv2[] */ + double tw; /* Total weight */ + int rc; + + if (smp[i].vflag == 0) + continue; + + /* Compute filtered value */ + tw = fdv2[0] = fdv2[1] = fdv2[2] = 0.0; + for (j = 0; j < smp[i].nnd; j++) { + nearsmth *np = smp[i].nd[j].n; /* Pointer to neighbor */ + double nw = smp[i].nd[j].w; /* Weight */ + double tmp[3]; + + if (np->vflag) { + icmSub3(tmp, smp[i].sv2, np->sv2); /* Vector from neighbour src to src */ + icmAdd3(tmp, tmp, np->dv2); /* Neigbour dst + vector */ + + icmScale3(tmp, tmp, nw); /* weight for filter */ + icmAdd3(fdv2, fdv2, tmp); /* sum filtered value */ + tw += nw; + } + } + + if (tw > 0.0) { +//printf("~1 %d: moved %f %f %f -> %f %f %f de %f\n", i, smp[i].dv2[0], smp[i].dv2[1], smp[i].dv2[2], fdv2[0], fdv2[1], fdv2[2], icmNorm33(smp[i].dv2,fdv2)); + icmScale3(smp[i].dv2, fdv2, 1.0/tw); + } + } +#endif /* SUBVEC_SMOOTHING */ + + VB(("near_smooth is done\n")); + + if (evectmap != NULL) + evectmap->del(evectmap); +#ifndef PLOT_DIGAM + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + for (i = 0; i < nmpts; i++) { + smp[i].sgam = NULL; + smp[i].dgam = NULL; + } +#else /* !PLOT_DIGAM */ + warning("!!!!! PLOT_DIGAM defined !!!!!"); +#endif /* !PLOT_DIGAM */ + *npp = nmpts; + return smp; +} + +/* Free the list of points that was returned */ +void free_nearsmth(nearsmth *smp, int nmpts) { + int i; + + for (i = 0; i < nmpts; i++) { + if (smp[i].nd != NULL) + free(smp[i].nd); + } + free(smp); +} + + +/* =================================================================== */ + +#if defined(SAVE_VRMLS) && defined(PLOT_MAPPING_INFLUENCE) +/* Create a plot indicating how the source mapping has been guided by the */ +/* various weighting forces */ +static void create_influence_plot(nearsmth *smp, int nmpts) { + gamut *gam; + int src = 0; /* 1 = src, 0 = dst gamuts */ + vrml *wrl = NULL; + co *fpnts = NULL; /* Mapping points to create diagnostic color mapping */ + rspl *swdiag = NULL; + int gres[3]; + double avgdev[3]; + double cols[4][3] = { { 1.0, 0.0, 0.0 }, /* Absolute = red */ + { 1.0, 1.0, 0.0 }, /* Relative = yellow */ + { 0.0, 0.0, 1.0 }, /* Radial = blue */ + { 0.0, 1.0, 0.0 } }; /* Depth = green */ + double grey[3] = { 0.5, 0.5, 0.5 }; /* Grey */ + double max, min; + int ix; + + if (src) + gam = sci_gam; + else + gam = di_gam; + + /* Setup the scattered data points */ + if ((fpnts = (co *)malloc((nmpts) * sizeof(co))) == NULL) { + fprintf(stderr,"gamut map: Malloc of diagnostic mapping setup points failed\n"); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + free_nearsmth(smp, nmpts); + *npp = 0; + return NULL; + } + + /* Compute error values and diagnostic color */ + /* for each guide vector */ + for (i = 0; i < nmpts; i++) { + double dv[4], gv; + double rgb[3]; + + /* Source value location */ + if (src) { + for (j = 0; j < 3; j++) + fpnts[i].p[j] = smp[i]._sv[j]; /* Non rotated and elevated */ + } else { /* Dest value location */ + for (j = 0; j < 3; j++) + fpnts[i].p[j] = smp[i].dv[j]; + } + + /* Diagnostic color */ + max = -1e60; min = 1e60; + for (k = 0; k < 4; k++) { /* Find max and min error value */ + dv[k] = smp[i].dbgv[k]; + if (dv[k] > max) + max = dv[k]; + if (dv[k] < min) + min = dv[k]; + } + for (k = 0; k < 4; k++) /* Scale to max */ + dv[k] /= max; + max /= max; + min /= max; + max -= min; /* reduce min to zero */ + for (k = 0; k < 4; k++) + dv[k] /= max; + for (gv = 1.0, k = 0; k < 4; k++) /* Blend remainder with grey */ + gv -= dv[k]; + + for (j = 0; j < 3; j++) /* Compute interpolated color */ + fpnts[i].v[j] = 0.0; + for (k = 0; k < 4; k++) { + for (j = 0; j < 3; j++) + fpnts[i].v[j] += dv[k] * cols[k][j]; + } + for (j = 0; j < 3; j++) + fpnts[i].v[j] += gv * grey[j]; + } + + /* Create the diagnostic color rspl */ + for (j = 0; j < 3; j++) { /* Set resolution for all axes */ + gres[j] = mapres; + avgdev[j] = 0.001; + } + swdiag = new_rspl(RSPL_NOFLAGS, 3, 3); /* Allocate 3D -> 3D */ + swdiag->fit_rspl(swdiag, RSPL_NOFLAGS, fpnts, nmpts, NULL, NULL, gres, NULL, NULL, 1.0, avgdev, NULL); + + /* Now create a plot of the sci_gam with the vertexes colored acording to the */ + /* diagnostic map. */ + if ((wrl = new_vrml("sci_gam_wt.wrl", 1)) == NULL) { + fprintf(stderr,"gamut map: new_vrml failed\n"); + if (fpnts != NULL) + free(fpnts); + if (swdiag != NULL) + swdiag->del(swdiag); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); + if (di_gam != sci_gam && di_gam != sci_gam) + di_gam->del(di_gam); + free_nearsmth(smp, nmpts); + *npp = 0; + return NULL; + } + + /* Plot the gamut triangle vertexes */ + for (ix = 0; ix >= 0;) { + co pp; + double col[3]; + + ix = gam->getvert(gam, NULL, pp.p, ix); + swdiag->interp(swdiag, &pp); + wrl->add_col_vertex(wrl, 0, pp.p, pp.v); + } + gam->startnexttri(gam); + for (;;) { + int vix[3]; + if (gam->getnexttri(gam, vix)) + break; + wrl->add_triangle(wrl, 0, vix); + } + wrl->make_triangles_vc(wrl, 0, 0.0); + + printf("Writing sci_gam_wt.wrl file\n"); + wrl->del(wrl); /* Write file */ + free(fpnts); + swdiag->del(swdiag); +} +#endif + + + + + + + + + + + + + + + + + + |