From 094535c010320967639e8e86f974d878e80baa72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Fri, 1 May 2015 16:13:57 +0200 Subject: Imported Upstream version 1.7.0 --- gamut/gammap.c | 325 ++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 220 insertions(+), 105 deletions(-) (limited to 'gamut/gammap.c') diff --git a/gamut/gammap.c b/gamut/gammap.c index d463755..638a2cd 100644 --- a/gamut/gammap.c +++ b/gamut/gammap.c @@ -45,20 +45,23 @@ #define VERBOSE /* [Def] Print out extra interesting information when verbose is set */ #undef PLOT_DIAG_WRL /* [Und] Always plot "gammap.wrl" */ - /* What do display when user requests disgnostic .wrl */ + /* What do display when user requests disgnostic VRML/X3D */ #define PLOT_SRC_GMT /* [Def] Plot the source surface to "gammap.wrl" as well */ #define PLOT_DST_GMT /* [Def] Plot the dest surface to "gammap.wrl" as well */ #undef PLOT_SRC_CUSPS /* [Und] Plot the source surface cusps to "gammap.wrl" as well */ #undef PLOT_DST_CUSPS /* [Und] Plot the dest surface cusps to "gammap.wrl" as well */ #undef PLOT_TRANSSRC_CUSPS /* [Und] Plot the gamut mapped source surface cusps to "gammap.wrl" */ -#undef PLOT_AXES /* [Und] Plot the axes to "gammap.wrl" as well */ +#define PLOT_AXES /* [Und] Plot the axes to "gammap.wrl" as well */ #undef SHOW_VECTOR_INDEXES /* [Und] Show the mapping vector index numbers */ -#define SHOW_MAP_VECTORS /* [Def] Show the mapping vectors */ -#undef SHOW_SUB_SURF /* [Und] Show the sub-surface mapping vector */ +#define SHOW_MAP_VECTORS /* [Def] Show the mapping vectors - yellow to red, green to red */ + /* if no clear direction. */ +#undef SHOW_SUB_SURF /* [Und] Show the sub-surface mapping vector - grey to purple. */ +#undef SHOW_SUB_PNTS /* [Und] Show the sub-surface sv2 (red), div2 (green), sd3 (yellow) pnts */ #undef SHOW_CUSPMAP /* [Und] Show the cusp mapped vectors rather than final vectors */ -#undef SHOW_ACTUAL_VECTORS /* [Und?] Show how the source vectors actually map thought xform */ +#undef SHOW_ACTUAL_VECTORS /* [Und] Show how the source vectors actually map thought xform */ #undef SHOW_ACTUAL_VEC_DIFF /* [Und] Show how the difference between guide and actual vectors */ + /* Other diagnostics */ #undef PLOT_LMAP /* [Und] Plot L map */ #undef PLOT_GAMUTS /* [Und] Save (part mapped) input and output gamuts as */ /* src.wrl, img.wrl, dst.wrl, gmsrc.wrl */ @@ -75,7 +78,11 @@ #undef PLOT_DIGAM /* [Und] Rather than DST_GMT - don't free it (#def in nearsmth.c too) */ -#define XRES 100 /* Res of plots */ +#define XRES 100 /* [100] Res of plots */ + +/* The locus.ts file can contain source locus(es) that will be plotted */ +/* as cones in red, with the destination plotted in white. They can */ +/* be created from .tif files using xicc/tiffgmts utility. */ /* Optional marker points for gamut mapping diagnosotic */ struct { @@ -85,6 +92,7 @@ struct { double col[3]; /* RGB color */ } markers[] = { { 0, }, /* End marker */ + { 1, { 37.18, 17.78, 20.28 }, { 0.545, 0.357, 0.256 } }, /* Dark Baby skin */ { 1, { 12.062, -0.87946, 0.97008 }, { 1.0, 0.3, 0.3 } }, /* Black point */ { 1, { 67.575411, -37.555250, -36.612862 }, { 1.0, 0.3, 0.3 } }, /* bad source in red (Red) */ { 1, { 61.003078, -44.466554, 1.922585 }, { 0.0, 1.0, 0.3 } }, /* good source in green */ @@ -169,7 +177,7 @@ static void inv_grey_func(void *pp, double *out, double *in); static void adjust_wb_func(void *pp, double *out, double *in); static void adjust_sat_func(void *pp, double *out, double *in); -#define XVRA 4.0 /* Extra mapping vertex ratio over no. tri verts from gamut */ +#define XVRA 3.0 /* [3.0] Extra mapping vertex ratio over no. tri verts from gamut */ /* The smoothed near weighting control values. */ /* These weightings setup the detailed behaviour of the */ @@ -186,16 +194,17 @@ gammapweights pweights[] = { 0.0, /* Cusp chroma alignment weighting 0 = none, 1 = full */ 0.3 /* Cusp hue alignment weighting 0 = none, 1 = full */ }, + 2.0, /* Alignment twist power, 0 = linear, 1 = curve, 2+ late curve */ 1.00 /* Chroma expansion 1 = none */ }, - { /* Radial weighting */ + { /* Radial weighting (currently broken - need to fix) */ 0.0, /* Radial error overall weight, 0 + */ 0.5, /* Radial hue dominance vs l+c, 0 - 1 */ 0.5 /* Radial l dominance vs, c, 0 - 1 */ }, { /* Weighting of absolute error of destination from source */ 1.0, /* Absolute error overall weight */ - 0.6, /* Hue dominance vs l+c, 0 - 1 */ + 0.8, /* Hue dominance vs l+c, 0 - 1 */ 0.8, /* White l dominance vs, c, 0 - 1 */ 0.5, /* Grey l dominance vs, c, 0 - 1 */ @@ -221,16 +230,16 @@ gammapweights pweights[] = { 0.0 /* Fine tuning expansion weight, 0 - 1 */ } }, -#ifdef NEVER { - gmm_l_d_blue, /* Increase maintaining hue importance for blue */ + gmm_light_yellow, /* Treat yellow differently, to get purer result. */ { { - -1.0, /* Cusp luminance alignment weighting 0 = none, 1 = full */ - -1.0, /* Cusp chroma alignment weighting 0 = none, 1 = full */ - 0.0 /* Cusp hue alignment weighting 0 = none, 1 = full */ + 0.9, /* Cusp luminance alignment weighting 0 = none, 1 = full */ + 0.8, /* Cusp chroma alignment weighting 0 = none, 1 = full */ + 0.7 /* Cusp hue alignment weighting 0 = none, 1 = full */ }, - -1.0 /* Chroma expansion 1 = none */ + 4.0, /* Alignment twist power, 0 = linear, 1 = curve, 2+ late curve */ + 1.20 /* Chroma expansion 1 = none */ }, { /* Radial weighting */ -1.0, /* Radial error overall weight, 0 + */ @@ -239,7 +248,7 @@ gammapweights pweights[] = { }, { /* Weighting of absolute error of destination from source */ -1.0, /* Absolute error overall weight */ - 0.8, /* Hue dominance vs l+c, 0 - 1 */ + -1.0, /* Hue dominance vs l+c, 0 - 1 */ -1.0, /* White l dominance vs, c, 0 - 1 */ -1.0, /* Grey l dominance vs, c, 0 - 1 */ @@ -252,7 +261,7 @@ gammapweights pweights[] = { -1.0 /* L error xover threshold in DE */ }, { /* Relative error preservation using smoothing */ - -1.0, 25.0 /* Relative Smoothing radius L* H* */ + 20.0, 10.0 /* Relative Smoothing radius L* H* */ }, { /* Weighting of excessive compression error, which is */ /* the src->dst vector length over the available dst depth. */ @@ -262,19 +271,20 @@ gammapweights pweights[] = { -1.0 /* Expansion depth weight */ }, { - -1.0 /* Fine tuning expansion weight, 0 - 1 */ + 0.5 /* Fine tuning expansion weight, 0 - 1 */ } }, -#endif /* NEVER */ +#ifdef NEVER { - gmm_light_yellow, /* Treat yellow differently, to get purer result. */ + gmm_l_d_blue, /* Increase maintaining hue importance for blue */ { { - 0.9, /* Cusp luminance alignment weighting 0 = none, 1 = full */ - 0.8, /* Cusp chroma alignment weighting 0 = none, 1 = full */ - 0.7 /* Cusp hue alignment weighting 0 = none, 1 = full */ + -1.0, /* Cusp luminance alignment weighting 0 = none, 1 = full */ + -1.0, /* Cusp chroma alignment weighting 0 = none, 1 = full */ + 0.0 /* Cusp hue alignment weighting 0 = none, 1 = full */ }, - 1.15 /* Chroma expansion 1 = none */ + -1.0, /* 2.0 Alignment twist power, 0 = linear, 1 = curve, 2+ late curve */ + -1.0 /* Chroma expansion 1 = none */ }, { /* Radial weighting */ -1.0, /* Radial error overall weight, 0 + */ @@ -283,7 +293,7 @@ gammapweights pweights[] = { }, { /* Weighting of absolute error of destination from source */ -1.0, /* Absolute error overall weight */ - -1.0, /* Hue dominance vs l+c, 0 - 1 */ + 0.8, /* Hue dominance vs l+c, 0 - 1 */ -1.0, /* White l dominance vs, c, 0 - 1 */ -1.0, /* Grey l dominance vs, c, 0 - 1 */ @@ -296,7 +306,7 @@ gammapweights pweights[] = { -1.0 /* L error xover threshold in DE */ }, { /* Relative error preservation using smoothing */ - 20.0, 20.0 /* Relative Smoothing radius L* H* */ + -1.0, 15.0 /* Relative Smoothing radius L* H* */ }, { /* Weighting of excessive compression error, which is */ /* the src->dst vector length over the available dst depth. */ @@ -306,14 +316,15 @@ gammapweights pweights[] = { -1.0 /* Expansion depth weight */ }, { - 0.5 /* Fine tuning expansion weight, 0 - 1 */ + -1.0 /* Fine tuning expansion weight, 0 - 1 */ } }, +#endif /* NEVER */ { gmm_end, } }; -double psmooth = 5.0; /* [5.0] Level of RSPL smoothing for perceptual, 1 = nominal */ +double psmooth = 4.0; /* [5.0] Level of RSPL smoothing for perceptual, 1 = nominal */ /* Saturation mapping weights, where saturation has priority over smoothness */ gammapweights sweights[] = { @@ -325,6 +336,7 @@ gammapweights sweights[] = { 0.5, /* Cusp chroma alignment weighting 0 = none, 1 = full */ 0.6 /* Cusp hue alignment weighting 0 = none, 1 = full */ }, + 1.0, /* Alignment twist power, 0 = linear, 1 = curve, 2+ late curve */ 1.05 /* Chroma expansion 1 = none */ }, { /* Radial weighting */ @@ -347,7 +359,7 @@ gammapweights sweights[] = { 10.0 /* L error extra xover threshold in DE */ }, { /* Relative vector smoothing */ - 15.0, 20.0 /* Relative Smoothing radius L* H* */ + 20.0, 25.0 /* Relative Smoothing radius L* H* */ }, { /* Weighting of excessive compression error, which is */ /* the src->dst vector length over the available dst depth. */ @@ -368,6 +380,7 @@ gammapweights sweights[] = { 1.0, /* Cusp chroma alignment weighting 0 = none, 1 = full */ 1.0 /* Cusp hue alignment weighting 0 = none, 1 = full */ }, + 1.0, /* Alignment twist power, 0 = linear, 1 = curve, 2+ late curve */ 1.20 /* Chroma expansion 1 = none */ }, { /* Radial weighting */ @@ -428,6 +441,7 @@ static void domap(gammap *s, double *out, double *in); static void dopartialmap1(gammap *s, double *out, double *in); static void dopartialmap2(gammap *s, double *out, double *in); static gamut *parttransgamut(gammap *s, gamut *src); +static void invdomap1(gammap *s, double *out, double *in); #ifdef PLOT_GAMUTS static void map_trans(void *cntx, double out[3], double in[3]); #endif @@ -487,7 +501,7 @@ gammap *new_gammap( int j; #if defined(PLOT_LMAP) || defined(PLOT_GAMUTS) || defined(PLOT_3DKNEES) - fprintf(stderr,"##### A gammap.c PLOT is #defined ####\n"); +# pragma message("################ A gammap.c PLOT is #defined #########################") #endif if (verb) { @@ -504,6 +518,7 @@ gammap *new_gammap( /* Setup methods */ s->del = del_gammap; s->domap = domap; + s->invdomap1 = invdomap1; /* Now create everything */ @@ -898,10 +913,50 @@ glumknf = 1.0; } } - /* To ensure symetry between compression and expansion, always create RSPL */ - /* for compression and its inverse, and then swap grey and igrey rspl to compensate. */ - if ((dwL - dbL) > (swL - sbL)) +#ifdef PLOT_LMAP + printf("sbL = %f, swL = %f\n",sbL,swL); + printf("dbL = %f, dwL = %f\n",dbL,dwL); +#endif + + /* Remember our source and destination mapping targets */ + /* so that we can use them for fine tuning later. */ + + /* We scale the source and target white and black */ + /* points to match the L values of the source and destination */ + /* L curve mapping, as this is how we have chosen the */ + /* white and black point mapping for the link. */ + /* Put them back in pre-rotated space, so that we can */ + /* check the overall transform of the white and black points. */ + t = (swL - sr_cs_bp[0])/(sr_cs_wp[0] - sr_cs_bp[0]); + for (j = 0; j < 3; j++) + s_mt_wp[j] = sr_cs_bp[j] + t * (sr_cs_wp[j] - sr_cs_bp[j]); + icmMul3By3x4(s_mt_wp, s->igrot, s_mt_wp); + + t = (sbL - sr_cs_wp[0])/(sr_cs_bp[0] - sr_cs_wp[0]); + for (j = 0; j < 3; j++) + s_mt_bp[j] = sr_cs_wp[j] + t * (sr_cs_bp[j] - sr_cs_wp[j]); +//printf("~1 check black point rotated = %f %f %f\n",s_mt_bp[0],s_mt_bp[1],s_mt_bp[2]); + icmMul3By3x4(s_mt_bp, s->igrot, s_mt_bp); +//printf("~1 check black point prerotated = %f %f %f\n",s_mt_bp[0],s_mt_bp[1],s_mt_bp[2]); + + t = (dwL - dr_cs_bp[0])/(dr_cs_wp[0] - dr_cs_bp[0]); + for (j = 0; j < 3; j++) + d_mt_wp[j] = dr_cs_bp[j] + t * (dr_cs_wp[j] - dr_cs_bp[j]); + + for (j = 0; j < 3; j++) + d_mt_bp[j] = dr_cs_wp[j] + t * (dr_cs_bp[j] - dr_cs_wp[j]); + + /* To ensure symetry between compression and expansion, always create RSPL for */ + /* overall compression and its inverse, and then swap grey and igrey rspl to compensate. */ + /* We swap the source and desitination white and black points to achieve this. */ + /* Note that we could still have expansion at one end or the other, depending */ + /* on the center point location, so we need to allow for this in the rspl setup. */ + if ((dwL - dbL) > (swL - sbL)) { + double tt; + tt = swL; swL = dwL; dwL = tt; + tt = sbL; sbL = dbL; dbL = tt; revrspl = 1; + } /* White point end */ lpnts[ngreyp].p[0] = swL; @@ -913,9 +968,6 @@ glumknf = 1.0; lpnts[ngreyp].v[0] = dbL; lpnts[ngreyp++].w = 10.0; /* Must go through here */ -//printf("~1 white loc %f, val %f\n",swL,dwL); -//printf("~1 black loc %f, val %f\n",sbL,dbL); - #ifdef USE_GLUMKNF if (gmi->glumknf < 0.05) #endif /* USE_GLUMKNF */ @@ -931,33 +983,42 @@ glumknf = 1.0; #ifdef USE_GLUMKNF else { /* There is at least some weight in knee points */ double cppos = 0.50; /* Center point ratio between black and white */ - double cplv; /* Center point location and value */ - double kppos = 0.30; /* Knee point ratio between white/black & center */ + double cpll, cplv; /* Center point location and value */ + double kpwpos = 0.30; /* White knee point location prop. towards center */ + double kpbpos = 0.20; /* Black knee point location prop. towards center */ double kwl, kbl, kwv, kbv; /* Knee point values and locations */ double kwx, kbx; /* Knee point extra */ +#ifdef PLOT_LMAP + printf("%ssbL = %f, swL = %f\n", revrspl ? "(swapped) ": "", sbL,swL); + printf("%sdbL = %f, dwL = %f\n", revrspl ? "(swapped) ": "", dbL,dwL); +#endif -//printf("sbL = %f, swL = %f\n",sbL,swL); -//printf("dbL = %f, dwL = %f\n",dbL,dwL); - - /* Center point */ + /* Center point location */ + cpll = cppos * (swL - sbL) + sbL; + // ~~?? would this be better if the output + // was scaled by dwL/swL ? cplv = cppos * (swL - sbL) + sbL; -//printf("~1 computed cplv = %f\n",cplv); +#ifdef PLOT_LMAP + printf("cpll = %f, cplv = %f\n",cpll, cplv); +#endif #ifdef NEVER /* Don't use a center point */ - lpnts[ngreyp].p[0] = cplv; + lpnts[ngreyp].p[0] = cpll; lpnts[ngreyp].v[0] = cplv; - lpnts[ngreyp++].w = 0.5; + lpnts[ngreyp++].w = 5.0; #endif //printf("~1 black half diff = %f\n",dbL - sbL); //printf("~1 white half diff = %f\n",dwL - swL); /* Knee point locations */ - kwl = kppos * (cplv - swL) + swL; - kbl = kppos * (cplv - sbL) + sbL; + kwl = kpwpos * (cplv - swL) + swL; + kbl = kpbpos * (cplv - sbL) + sbL; /* Extra compression for white and black knees */ + // ~~ ie move knee point level beyond 45 degree line + // ~~ weigting of black point and white point differences kwx = 0.6 * (dbL - sbL) + 1.0 * (swL - dwL); kbx = 1.0 * (dbL - sbL) + 0.6 * (swL - dwL); @@ -975,47 +1036,22 @@ glumknf = 1.0; kbv = dbL; -//printf("~1 kbl = %f, kbv = %f\n",kbl, kbv); -//printf("~1 kwl = %f, kwv = %f\n",kwl, kwv); +#ifdef PLOT_LMAP + printf("using kbl = %f, kbv = %f\n",kbl, kbv); + printf("using kwl = %f, kwv = %f\n",kwl, kwv); +#endif - /* Emphasise points to cause "knee" curve */ + /* Emphasise points to cause white "knee" curve */ lpnts[ngreyp].p[0] = kwl; lpnts[ngreyp].v[0] = kwv; lpnts[ngreyp++].w = gmi->glumknf * gmi->glumknf; + /* Emphasise points to cause black "knee" curve */ lpnts[ngreyp].p[0] = kbl; lpnts[ngreyp].v[0] = kbv; lpnts[ngreyp++].w = 1.5 * gmi->glumknf * 1.5 * gmi->glumknf; } #endif /* USE_GLUMKNF */ - - /* Remember our source and destinatio mapping targets */ - /* so that we can use them for fine tuning later. */ - - /* We scale the source and target white and black */ - /* points to match the L values of the source and destination */ - /* L curve mapping, as this is how we have chosen the */ - /* white and black point mapping for the link. */ - /* Put them back in pre-rotated space, so that we can */ - /* check the overall transform of the white and black points. */ - t = (swL - sr_cs_bp[0])/(sr_cs_wp[0] - sr_cs_bp[0]); - for (j = 0; j < 3; j++) - s_mt_wp[j] = sr_cs_bp[j] + t * (sr_cs_wp[j] - sr_cs_bp[j]); - icmMul3By3x4(s_mt_wp, s->igrot, s_mt_wp); - - t = (sbL - sr_cs_wp[0])/(sr_cs_bp[0] - sr_cs_wp[0]); - for (j = 0; j < 3; j++) - s_mt_bp[j] = sr_cs_wp[j] + t * (sr_cs_bp[j] - sr_cs_wp[j]); -//printf("~1 check black point rotated = %f %f %f\n",s_mt_bp[0],s_mt_bp[1],s_mt_bp[2]); - icmMul3By3x4(s_mt_bp, s->igrot, s_mt_bp); -//printf("~1 check black point prerotated = %f %f %f\n",s_mt_bp[0],s_mt_bp[1],s_mt_bp[2]); - - t = (dwL - dr_cs_bp[0])/(dr_cs_wp[0] - dr_cs_bp[0]); - for (j = 0; j < 3; j++) - d_mt_wp[j] = dr_cs_bp[j] + t * (dr_cs_wp[j] - dr_cs_bp[j]); - - for (j = 0; j < 3; j++) - d_mt_bp[j] = dr_cs_wp[j] + t * (dr_cs_bp[j] - dr_cs_wp[j]); } /* We now create the 1D rspl L map, that compresses or expands the luminence */ @@ -1031,15 +1067,6 @@ glumknf = 1.0; double avgdev[MXDO]; int gres = 256; - if (revrspl) { /* Invert creation and usage for symetry between compress and exp. */ - int i; - for (i = 0; i < ngreyp; i++) { - double tt = lpnts[i].p[0]; /* Swap source and dest */ - lpnts[i].p[0] = lpnts[i].v[0]; - lpnts[i].v[0] = tt; - } - } - /* Create a 1D rspl, that is used to */ /* form the overall L compression mapping. */ if ((s->grey = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) /* Allocate 1D -> 1D */ @@ -1076,7 +1103,7 @@ glumknf = 1.0; /* Create it from inverse lookups of s->grey */ s->igrey->set_rspl(s->igrey, 0, (void *)s->grey, inv_grey_func, il, ih, &gres, ol, oh); - if (revrspl) { /* Swap to compensate for expansion */ + if (revrspl) { /* Swap to compensate for swapping of white and black points */ rspl *tt = s->grey; s->grey = s->igrey; s->igrey = tt; @@ -1140,6 +1167,7 @@ glumknf = 1.0; /* if there is any compression or expansion to do. */ if (gmi->gamcpf > 1e-6 || gmi->gamexf > 1e-6) { cow *gpnts = NULL; /* Mapping points to create gamut mapping */ + int max_gpnts; int nspts; /* Number of source gamut surface points */ int rgridpts; /* Number of range surface grid points */ int i, j; @@ -1175,7 +1203,8 @@ typedef struct { } #endif - if ((gpnts = (cow *)malloc((nres + 3 * nspts + rgridpts) * sizeof(cow))) == NULL) { + max_gpnts = nres + 3 * nspts + rgridpts; + if ((gpnts = (cow *)malloc(max_gpnts * sizeof(cow))) == NULL) { fprintf(stderr,"gamut map: Malloc of mapping setup points failed\n"); s->grey->del(s->grey); s->igrey->del(s->igrey); @@ -1275,6 +1304,8 @@ typedef struct { gpnts[ngamp].w); #endif ngamp++; + if (ngamp >= max_gpnts) + error("gammap: internal, not enough space for mapping points A (%d > %d)\n",ngamp, max_gpnts); } #endif /* USE_GREYMAP */ @@ -1436,10 +1467,10 @@ typedef struct { doaxes = 1; #endif if (diagname != NULL) - wrl = new_vrml(diagname, doaxes, 0); + wrl = new_vrml(diagname, doaxes, vrml_lab); #ifdef PLOT_DIAG_WRL else - wrl = new_vrml("gammap.wrl", doaxes, 0); + wrl = new_vrml("gammap", doaxes, vrml_lab); #endif } @@ -1513,31 +1544,51 @@ typedef struct { } gpnts[ngamp++].w = 1.01; /* Main gamut surface mapping point */ /* (Use 1.01 as a marker value) */ + if (ngamp >= max_gpnts) + error("gammap: internal, not enough space for mapping points B (%d > %d)\n",ngamp, max_gpnts); #ifdef USE_GAMKNF /* Add sub surface mapping point if available */ if (nsm[i].vflag != 0) { /* Sub surface point is available */ /* Compute destination value which is a blend */ - /* between the source value and the fully mapped destination value. */ + /* between the source value and the knee adjusted destination */ icmBlend3(nsm[i].div2, nsm[i].sv2, nsm[i].dv2, cpexf); #ifdef NEVER printf("Src2 point = %f %f %f radius %f\n",nsm[i].sv2[0], nsm[i].sv2[1], nsm[i].sv2[2], nsm[i].sr); printf("Dst2 point = %f %f %f radius %f\n",nsm[i].dv2[0], nsm[i].dv2[1], nsm[i].dv2[2], nsm[i].dr); printf("Blended dst2 point = %f %f %f\n",nsm[i].div2[0], nsm[i].div2[1], nsm[i].div2[2]); + printf("Src/Dst3 point = %f %f %f w %f\n",nsm[i].sd2[0], nsm[i].sd2[1], nsm[i].sd2[2]); printf("\n"); #endif /* NEVER */ + /* Set the sub-surface gamut hull mapping point */ for (j = 0; j < 3; j++) { gpnts[ngamp].p[j] = nsm[i].sv2[j]; gpnts[ngamp].v[j] = nsm[i].div2[j]; } gpnts[ngamp++].w = nsm[i].w2; /* Sub-suface mapping points */ + + if (ngamp >= max_gpnts) + error("gammap: internal, not enough space for mapping points C (%d > %d)\n",ngamp, max_gpnts); + + /* Set the sub-surface gamut hull mapping point */ + for (j = 0; j < 3; j++) { + gpnts[ngamp].p[j] = nsm[i].sd3[j]; + gpnts[ngamp].v[j] = nsm[i].sd3[j]; + } + gpnts[ngamp++].w = nsm[i].w3; /* Sub-suface mapping points */ + + if (ngamp >= max_gpnts) + error("gammap: internal, not enough space for mapping points D (%d > %d)\n",ngamp, max_gpnts); } #endif /* USE_GAMKNF */ } + if (ngamp >= max_gpnts) + error("gammap: internal, not enough space for mapping points (%d > %d)\n",ngamp, max_gpnts); + /* Create preliminary gamut mapping rspl, without grid boundary values. */ /* We use this to lookup the mapping for points on the source space gamut */ /* that result from clipping our grid boundary points */ @@ -1601,6 +1652,9 @@ typedef struct { /* A low weight seems to be enough ? */ /* the lower the better in terms of geting best hull mapping fidelity */ gpnts[ngamp++].w = 0.05 * ww; + + if (ngamp >= max_gpnts) + error("gammap: internal, not enough space for mapping points E (%d > %d)\n",ngamp, max_gpnts); } DC_INC(gc); if (DC_DONE(gc)) @@ -1833,11 +1887,11 @@ typedef struct { /* Add the source and dest gamut surfaces */ #ifdef PLOT_SRC_GMT - wrl->make_gamut_surface_2(wrl, sil_gam, 0.6, 0, cc); + wrl->make_gamut_surface_2(wrl, sil_gam, 0.6, 0, cc); /* Grey */ #endif /* PLOT_SRC_GMT */ #ifdef PLOT_DST_GMT cc[0] = -1.0; - wrl->make_gamut_surface(wrl, d_gam, 0.2, cc); + wrl->make_gamut_surface(wrl, d_gam, 0.3, cc); /* Natural color */ #endif /* PLOT_DST_GMT */ #ifdef PLOT_DIGAM if (nsm[0].dgam == NULL) @@ -1918,7 +1972,6 @@ typedef struct { # ifdef SHOW_SUB_SURF if (nsm[i].vflag != 0) { /* Sub surface point is available */ - wrl->add_col_vertex(wrl, 0, nsm[i].sv2, lgrey); /* Subs-surf Source value */ wrl->add_col_vertex(wrl, 0, nsm[i].div2, purp); /* Blended destination value */ } @@ -1928,14 +1981,29 @@ typedef struct { wrl->make_lines(wrl, 0, 2); /* Guide vectors */ #endif /* Show vectors */ -#ifdef SHOW_VECTOR_INDEXES +#if defined(SHOW_VECTOR_INDEXES) || defined(SHOW_SUB_PNTS) for (i = 0; i < nnsm; i++) { - double cream[3] = { 0.7, 0.7, 0.5 }; - char buf[100]; - sprintf(buf, "%d", i); - wrl->add_text(wrl, buf, nsm[i].sv, cream, 0.5); - } +#ifdef SHOW_VECTOR_INDEXES + { + double cream[3] = { 0.7, 0.7, 0.5 }; + char buf[100]; + sprintf(buf, "%d", i); + wrl->add_text(wrl, buf, nsm[i].sv, cream, 0.5); + } #endif /* SHOW_VECTOR_INDEXES */ +# ifdef SHOW_SUB_PNTS + if (nsm[i].vflag != 0) { /* Sub surface point is available */ + double red[3] = { 1.0, 0.0, 0.0 }; + double green[3] = { 0.0, 1.0, 0.0 }; + double yellow[3] = { 1.0, 1.0, 0.0 }; + + wrl->add_marker(wrl, nsm[i].sv2, red, 1.0); /* Subs-surf Source value */ + wrl->add_marker(wrl, nsm[i].div2, green, 1.0); /* Blended destination value */ + wrl->add_marker(wrl, nsm[i].sd3, yellow, 1.0); /* Deep sub-surface point */ + } +# endif /* SHOW_SUB_PNTS */ + } +#endif /* add the locus from locus.ts file */ if (locus != NULL) { @@ -2184,10 +2252,10 @@ typedef struct { } #ifdef PLOT_GAMUTS - scl_gam->write_vrml(scl_gam, "src.wrl", 1, 0); - sil_gam->write_vrml(sil_gam, "img.wrl", 1, 0); - d_gam->write_vrml(d_gam, "dst.wrl", 1, 0); - sc_gam->write_trans_vrml(sc_gam, "gmsrc.wrl", 1, 0, map_trans, s); + scl_gam->write_vrml(scl_gam, "src", 1, 0); + sil_gam->write_vrml(sil_gam, "img", 1, 0); + d_gam->write_vrml(d_gam, "dst", 1, 0); + sc_gam->write_trans_vrml(sc_gam, "gmsrc", 1, 0, map_trans, s); #endif if (sil_gam != scl_gam) @@ -2222,6 +2290,8 @@ gammap *s free(s); } +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + /* Apply the gamut mapping to the given color value */ static void domap( gammap *s, @@ -2336,6 +2406,51 @@ double *in } } +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/* powell function - minimise error to target */ +static double invgmfunc( +void *fdata, +double *tp +) { + gammap *s = (gammap *)fdata; + int i; + double gmv[3]; + double tt, rv = 0.0; + + domap(s, gmv, tp); + for (i = 0; i < 3; i++) { + double tt = gmv[i] - s->tv[i]; + rv += tt * tt; + } + + return rv; +} + +/* Invert a gamut mapping using powell */ +static void invdomap1( +gammap *s, +double *out, +double *in +) { + 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]; + + if (powell(&rv, 3, tp, ss, 1e-7, 5000, invgmfunc, (void *)s, NULL, NULL) != 0) { + warning("gamut invdomap1 failed on %f %f %f\n", in[0], in[1], in[2]); + } + + out[0] = tp[0]; + out[1] = tp[1]; + out[2] = tp[2]; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + /* Function to pass to rspl to invert grey curve */ static void inv_grey_func( void *cntx, -- cgit v1.2.3