diff options
Diffstat (limited to 'gamut')
| -rw-r--r-- | gamut/gammap.c | 667 | ||||
| -rw-r--r-- | gamut/gamut.c | 682 | ||||
| -rw-r--r-- | gamut/gamut.h | 25 | ||||
| -rw-r--r-- | gamut/maptest.c | 3 | ||||
| -rw-r--r-- | gamut/nearsmth.c | 1803 | ||||
| -rw-r--r-- | gamut/nearsmth.h | 48 | ||||
| -rw-r--r-- | gamut/smthtest.c | 35 |
7 files changed, 1358 insertions, 1905 deletions
diff --git a/gamut/gammap.c b/gamut/gammap.c index e276b99..638a2cd 100644 --- a/gamut/gammap.c +++ b/gamut/gammap.c @@ -65,9 +65,13 @@ #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 */ -#undef PLOT_3DKNEES /* [Und] Plot each 3D compression knee */ +#undef PLOT_3DKNEES /* [Und] Plot the 3D compression knees */ #undef CHECK_NEARMAP /* [Und] Check how accurately near map vectors are represented by rspl */ -#undef DUMP_GREY_AXIS_POINTS /* [Und] Dump grey axis map 3d->3d points */ + +#define USE_GLUMKNF /* [Define] Enable luminence knee function points */ +#define USE_GREYMAP /* [Define] Enable 3D->3D mapping points down the grey axis */ +#define USE_GAMKNF /* [Define] Enable 3D knee function points */ +#define USE_BOUND /* [Define] Enable grid boundary anchor points */ #undef SHOW_NEIGBORS /* [Und] Show nearsmth neigbors in gammap.wrl */ @@ -76,13 +80,6 @@ #define XRES 100 /* [100] Res of plots */ - /* Functionality */ -#define USE_GLUMKNF /* [Define] Enable luminence knee function points else linear */ -#define USE_GREYMAP /* [Define] Enable 3D->3D mapping points down the grey axis */ -#define USE_GAMKNF /* [Define] Enable 3D knee function points */ -#define USE_BOUND /* [Define] Enable grid boundary anchor points */ - - /* 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. */ @@ -151,6 +148,7 @@ struct { #include <fcntl.h> #include <string.h> #include <math.h> +#include "counters.h" #include "icc.h" #include "numlib.h" #include "xicc.h" @@ -170,17 +168,11 @@ typedef struct { double satenh; /* Saturation engancement value */ } adjustsat; -/* Callback context for fixing 1D Lut white and black points */ -typedef struct { - double twb[2], awb[2]; /* Target and Actual black, white */ -} adjust1wb; - /* Callback context for making clut relative to white and black points */ typedef struct { double mat[3][4]; } adjustwb; -static void adjust1_wb_func(void *pp, double *out, double *in); 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); @@ -200,7 +192,7 @@ gammapweights pweights[] = { { 0.1, /* Cusp luminance alignment weighting 0 = none, 1 = full */ 0.0, /* Cusp chroma alignment weighting 0 = none, 1 = full */ - 0.2 /* Cusp hue 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 */ @@ -214,9 +206,9 @@ gammapweights pweights[] = { 1.0, /* Absolute error overall weight */ 0.8, /* Hue dominance vs l+c, 0 - 1 */ - 0.8, /* [0.8] White l dominance vs, c, 0 - 1 */ - 0.45, /* [0.45] Grey l dominance vs, c, 0 - 1 */ - 0.94, /* [0.94] Black l dominance vs, c, 0 - 1 */ + 0.8, /* White l dominance vs, c, 0 - 1 */ + 0.5, /* Grey l dominance vs, c, 0 - 1 */ + 0.97, /* Black l dominance vs, c, 0 - 1 */ 0.4, /* White l blend start radius, 0 - 1, at white = 0 */ 0.7, /* Black l blend power, linear = 1.0, enhance < 1.0 */ @@ -225,15 +217,14 @@ gammapweights pweights[] = { 10.0 /* L error extra xover threshold in DE */ }, { /* Relative vector smoothing */ - 20.0, 30.0, /* Relative Smoothing radius L* H* */ - 0.9 /* Degree of smoothing */ + 25.0, 35.0 /* Relative Smoothing radius L* H* */ }, { /* Weighting of excessive compression error, which is */ /* the src->dst vector length over the available dst depth. */ /* The depth is half the distance to the intersection of the */ /* vector to the other side of the gamut. (doesn't get triggered much ?) */ - 5.0, /* [5] Compression depth weight */ - 5.0 /* [5] Expansion depth weight */ + 10.0, /* Compression depth weight */ + 10.0 /* Expansion depth weight */ }, { 0.0 /* Fine tuning expansion weight, 0 - 1 */ @@ -270,8 +261,7 @@ gammapweights pweights[] = { -1.0 /* L error xover threshold in DE */ }, { /* Relative error preservation using smoothing */ - 20.0, 10.0, /* Relative Smoothing radius L* H* */ - 0.5 /* Degree of smoothing */ + 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. */ @@ -303,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 */ @@ -316,8 +306,7 @@ gammapweights pweights[] = { -1.0 /* L error xover threshold in DE */ }, { /* Relative error preservation using smoothing */ - -1.0, 15.0, /* Relative Smoothing radius L* H* */ - -1.0 /* Degree of smoothing */ + -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. */ @@ -335,61 +324,7 @@ gammapweights pweights[] = { gmm_end, } }; -double psmooth = 2.0; /* [2.0] Level of RSPL smoothing for perceptual, 1 = nominal */ - -/* Lightness Preserving Perceptual mapping weights, where preserving lightness */ -/* and hue has the highest priority, followed by smoothness and proportionality. */ -/* Chroma is basically sacrificed. */ -gammapweights lpweights[] = { - { - gmm_default, /* Non hue specific defaults */ - { /* Cusp alignment control */ - { - 0.0, /* [0.2] Cusp luminance alignment weighting 0 = none, 1 = full */ - 0.0, /* [0.0] Cusp chroma alignment weighting 0 = none, 1 = full */ - 0.0 /* [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 (currently broken - need to fix) */ - 0.0, /* Radial error overall weight, 0 + */ - 0.1, /* Radial hue dominance vs l+c, 0 - 1 */ - 1.0 /* Radial l dominance vs, c, 0 - 1 */ - }, - { /* Weighting of absolute error of destination from source */ - 1.0, /* Absolute error overall weight */ - 0.90, /* [0.9] Hue dominance vs l+c, 0 - 1 */ - - 0.98, /* White l dominance vs, c, 0 - 1 */ - 0.95, /* Grey l dominance vs, c, 0 - 1 */ - 0.99, /* Black l dominance vs, c, 0 - 1 */ - - 0.4, /* White l blend start radius, 0 - 1, at white = 0 */ - 0.7, /* Black l blend power, linear = 1.0, enhance < 1.0 */ - - 1.0, /* L error extra power with size, none = 1.0 */ - 100.0 /* L error extra xover threshold in DE */ - }, - { /* Relative vector smoothing */ - 6.0, 30.0, /* Relative Smoothing radius L* H* */ - 0.9 /* [0.9] Degree of smoothing */ - }, - { /* Weighting of excessive compression error, which is */ - /* the src->dst vector length over the available dst depth. */ - /* (This compromizes constanl L near white and black, so minimize) */ - 0.0, /* Compression depth weight */ - 0.0 /* Expansion depth weight */ - }, - { - 0.0 /* Fine tuning expansion weight, 0 - 1 */ - } - }, - { - gmm_end, - } -}; -double lpsmooth = 1.0; /* [1.0] Level of RSPL smoothing for ligtness pres perc, 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[] = { @@ -414,25 +349,24 @@ gammapweights sweights[] = { 0.4, /* Hue dominance vs l+c, 0 - 1 */ 0.6, /* White l dominance vs, c, 0 - 1 */ - 0.3, /* Grey l dominance vs, c, 0 - 1 */ - 0.7, /* Black l dominance vs, c, 0 - 1 */ + 0.4, /* Grey l dominance vs, c, 0 - 1 */ + 0.6, /* Black l dominance vs, c, 0 - 1 */ 0.5, /* wl blend start radius, 0 - 1 */ 1.0, /* bl blend power, linear = 1.0, enhance < 1.0 */ - 1.5, /* L error extra power with size, none = 1.0 */ - 20.0 /* L error extra xover threshold in DE */ + 1.0, /* L error extra power with size, none = 1.0 */ + 10.0 /* L error extra xover threshold in DE */ }, { /* Relative vector smoothing */ - 15.0, 20.0, /* Relative Smoothing radius L* H* */ - 0.8 /* [0.8] Degree of smoothing */ + 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. */ /* The depth is half the distance to the intersection of the */ /* vector to the other side of the gamut. (doesn't get triggered much ?) */ - 5.0, /* Compression depth weight */ - 5.0 /* Expansion depth weight */ + 10.0, /* Compression depth weight */ + 10.0 /* Expansion depth weight */ }, { 0.5 /* Fine tuning expansion weight, 0 - 1 */ @@ -469,8 +403,7 @@ gammapweights sweights[] = { -1.0 /* L error xover threshold in DE */ }, { /* Relative error preservation using smoothing */ - 10.0, 15.0, /* Relative smoothing radius */ - 0.5 /* Degree of smoothing */ + 10.0, 15.0 /* Relative smoothing radius */ }, { /* Weighting of excessive compression error, which is */ /* the src->dst vector length over the available dst depth. */ @@ -487,9 +420,7 @@ gammapweights sweights[] = { gmm_end } }; -/* The cusp alignment tends to upset the vector smoothing (not exactly sure why), */ -/* so use more rspl smoothing to compensate. */ -double ssmooth = 4.0; /* [1.0] Level of RSPL smoothing for saturation */ +double ssmooth = 2.0; /* Level of RSPL smoothing for saturation */ /* * Notes: @@ -520,7 +451,7 @@ static void map_trans(void *cntx, double out[3], double in[3]); gammap *new_gammap( int verb, /* Verbose flag */ gamut *sc_gam, /* Source colorspace gamut */ - gamut *isi_gam, /* Input source image gamut (NULL if none) */ + gamut *si_gam, /* Source image gamut (NULL if none) */ gamut *d_gam, /* Destination colorspace gamut */ icxGMappingIntent *gmi, /* Gamut mapping specification */ int src_kbp, /* Use K only black point as src gamut black point */ @@ -532,10 +463,9 @@ gammap *new_gammap( double *mx, /* for rspl grid. */ char *diagname /* If non-NULL, write a gamut mapping diagnostic WRL */ ) { - gammap *s; /* This */ - gamut *si_gam = NULL; /* Source image gamut (intersected with sc_gam) */ - gamut *scl_gam; /* Source colorspace gamut with rotation and L mapping applied */ - gamut *sil_gam; /* Source image gamut with rotation and L mapping applied */ + gammap *s; /* This */ + gamut *scl_gam; /* Source colorspace gamut with rotation and L mapping applied */ + gamut *sil_gam; /* Source image gamut with rotation and L mapping applied */ double s_cs_wp[3]; /* Source colorspace white point */ double s_cs_bp[3]; /* Source colorspace black point */ @@ -561,11 +491,7 @@ gammap *new_gammap( double d_mt_wp[3]; /* Overall destination mapping white point (used for finetune) */ double d_mt_bp[3]; /* Overall destination mapping black point (used for finetune) */ -#ifdef USE_BOUND - int surfpnts = 1; /* Add grid surface anchor points */ -#else - int surfpnts = 0; /* Don't add grid surface points */ -#endif + int defrgrid = 6; /* mapping range surface default anchor point resolution */ int nres = 512; /* Neutral axis resolution */ cow lpnts[10]; /* Mapping points to create grey axis map */ int revrspl = 0; /* Reverse grey axis rspl construction */ @@ -578,14 +504,10 @@ gammap *new_gammap( # pragma message("################ A gammap.c PLOT is #defined #########################") #endif -#ifndef USE_BOUND -# pragma message("################ gammap.c USE_BOUND not set #########################") -#endif - if (verb) { xicc_dump_gmi(gmi); printf("Gamut map resolution: %d\n",mapres); - if (isi_gam != NULL) + if (si_gam != NULL) printf("Image gamut supplied\n"); } @@ -600,15 +522,18 @@ gammap *new_gammap( /* Now create everything */ - /* Grab the colorspace white and black points */ + /* Grab the white and black points */ if (src_kbp) { - if (sc_gam->getwb(sc_gam, s_cs_wp, NULL, s_cs_bp, NULL, NULL, NULL)) { + // ~~99 Hmm. Shouldn't this be colorspace rather than gamut ???? + if (sc_gam->getwb(sc_gam, NULL, NULL, NULL, s_cs_wp, NULL, s_cs_bp)) { +// if (sc_gam->getwb(sc_gam, s_cs_wp, NULL, s_cs_bp, NULL, NULL, NULL)) fprintf(stderr,"gamut map: Unable to read source colorspace white and black points\n"); free(s); return NULL; } } else { - if (sc_gam->getwb(sc_gam, s_cs_wp, s_cs_bp, NULL, NULL, NULL, NULL)) { + if (sc_gam->getwb(sc_gam, NULL, NULL, NULL, s_cs_wp, s_cs_bp, NULL)) { +// if (sc_gam->getwb(sc_gam, s_cs_wp, s_cs_bp, NULL, NULL, NULL, NULL)) fprintf(stderr,"gamut map: Unable to read source colorspace white and black points\n"); free(s); return NULL; @@ -616,40 +541,16 @@ gammap *new_gammap( } /* If source space is source gamut */ - if (isi_gam == NULL || isi_gam == sc_gam) { + if (si_gam == NULL) { si_gam = sc_gam; for (j = 0; j < 3; j++) { s_ga_wp[j] = s_cs_wp[j]; s_ga_bp[j] = s_cs_bp[j]; } - /* Else have explicit image gamut */ + /* Else have explicit sourcegamut */ } else { -#ifdef VERBOSE - if (verb) { /* Check that image gamut is within colorspace */ - double scwp[3], scbp[3]; - double imwp[3], imbp[3]; - - sc_gam->getwb(sc_gam, NULL, NULL, NULL, scwp, scbp, NULL); - isi_gam->getwb(isi_gam, NULL, NULL, NULL, imwp, imbp, NULL); - - if (imwp[0] > (scwp[0] + 1e-4) - || imbp[0] < (scbp[0] - 1e-4)) { - printf("Warning: image gamut is bigger than src colorspace!\n"); - } - - } -#endif - /* Intersect it with the source colorspace gamut in case */ - /* something strange is going on. (mismatched appearance params ?) */ - if ((si_gam = new_gamut(0.0, 0, 0)) == NULL) { - fprintf(stderr,"gamut map: new_gamut failed\n"); - free(s); - return NULL; - } - si_gam->intersect(si_gam, isi_gam, sc_gam); - if (src_kbp) { if (si_gam->getwb(si_gam, NULL, NULL, NULL, s_ga_wp, NULL, s_ga_bp)) { fprintf(stderr,"gamut map: Unable to read source gamut white and black points\n"); @@ -663,21 +564,42 @@ gammap *new_gammap( return NULL; } } + + /* Guard against silliness. Image must be within colorspace */ + if (s_ga_wp[0] > s_cs_wp[0]) { + int j; + double t; +#ifdef VERBOSE + if (verb) + printf("Fixing wayward image white point\n"); +#endif + t = (s_cs_wp[0] - s_ga_bp[0])/(s_ga_wp[0] - s_ga_bp[0]); + for (j = 0; j < 3; j++) + s_ga_wp[j] = s_ga_bp[j] + t * (s_ga_wp[j] - s_ga_bp[j]); + + } + if (s_ga_bp[0] < s_cs_bp[0]) { + int j; + double t; +#ifdef VERBOSE + if (verb) + printf("Fixing wayward image black point\n"); +#endif + t = (s_cs_bp[0] - s_ga_wp[0])/(s_ga_bp[0] - s_ga_wp[0]); + for (j = 0; j < 3; j++) + s_ga_bp[j] = s_ga_wp[j] + t * (s_ga_bp[j] - s_ga_wp[j]); + } } if (dst_kbp) { if (d_gam->getwb(d_gam, NULL, NULL, NULL, d_cs_wp, NULL, d_cs_bp)) { fprintf(stderr,"gamut map: Unable to read destination white and black points\n"); - if (si_gam != sc_gam) - si_gam->del(si_gam); free(s); return NULL; } } else { if (d_gam->getwb(d_gam, NULL, NULL, NULL, d_cs_wp, d_cs_bp, NULL)) { fprintf(stderr,"gamut map: Unable to read destination white and black points\n"); - if (si_gam != sc_gam) - si_gam->del(si_gam); free(s); return NULL; } @@ -950,7 +872,7 @@ glumknf = 1.0; /* If we have a gamut (ie. image) range that is smaller than the */ /* L range of the colorspace, then use its white and black L values */ /* as the source to be compressed to the destination L range. */ - /* We expand only a colorspace range, not an image gamut range. */ + /* We expand only a colorspace range, not a gamut/image range. */ { double swL, dwL; /* Source and destination white point L */ double sbL, dbL; /* Source and destination black point L */ @@ -963,14 +885,14 @@ glumknf = 1.0; dwL = gmi->glumwexf * dr_cs_wp[0] + (1.0 - gmi->glumwexf) * sr_cs_wp[0]; } else { - if (sr_cs_wp[0] > dr_cs_wp[0]) { /* Colorspace needs compression */ + if (sr_ga_wp[0] > dr_cs_wp[0]) { /* Gamut or colorspace needs compression */ - swL = (1.0 - gmi->glumwcpf) * dr_cs_wp[0] + gmi->glumwcpf * sr_cs_wp[0]; + swL = (1.0 - gmi->glumwcpf) * dr_cs_wp[0] + gmi->glumwcpf * sr_ga_wp[0]; dwL = dr_cs_wp[0]; } else { /* Neither needed */ - swL = sr_cs_wp[0]; - dwL = sr_cs_wp[0]; + swL = sr_ga_wp[0]; + dwL = sr_ga_wp[0]; } } @@ -980,14 +902,14 @@ glumknf = 1.0; dbL = gmi->glumbexf * dr_cs_bp[0] + (1.0 - gmi->glumbexf) * sr_cs_bp[0]; } else { - if (sr_cs_bp[0] < dr_cs_bp[0]) { /* Colorspace needs compression */ + if (sr_ga_bp[0] < dr_cs_bp[0]) { /* Gamut or colorspace needs compression */ - sbL = (1.0 - gmi->glumbcpf) * dr_cs_bp[0] + gmi->glumbcpf * sr_cs_bp[0]; + sbL = (1.0 - gmi->glumbcpf) * dr_cs_bp[0] + gmi->glumbcpf * sr_ga_bp[0]; dbL = dr_cs_bp[0]; } else { /* Neither needed */ - sbL = sr_cs_bp[0]; - dbL = sr_cs_bp[0]; + sbL = sr_ga_bp[0]; + dbL = sr_ga_bp[0]; } } @@ -1046,21 +968,24 @@ glumknf = 1.0; lpnts[ngreyp].v[0] = dbL; lpnts[ngreyp++].w = 10.0; /* Must go through here */ -#ifndef USE_GLUMKNF - /* make sure curve is firmly anchored */ - lpnts[ngreyp].p[0] = 0.3 * lpnts[ngreyp-1].p[0] + 0.7 * lpnts[ngreyp-2].p[0]; - lpnts[ngreyp].v[0] = 0.3 * lpnts[ngreyp-1].v[0] + 0.7 * lpnts[ngreyp-2].v[0]; - lpnts[ngreyp++].w = 1.0; - - lpnts[ngreyp].p[0] = 0.7 * lpnts[ngreyp-2].p[0] + 0.3 * lpnts[ngreyp-3].p[0]; - lpnts[ngreyp].v[0] = 0.7 * lpnts[ngreyp-2].v[0] + 0.3 * lpnts[ngreyp-3].v[0]; - lpnts[ngreyp++].w = 1.0; -#else /* USE_GLUMKNF */ - { - double cppos = 0.50; /* [0.50] Center point ratio between black and white */ +#ifdef USE_GLUMKNF + if (gmi->glumknf < 0.05) +#endif /* USE_GLUMKNF */ + { /* make sure curve is firmly anchored */ + lpnts[ngreyp].p[0] = 0.3 * lpnts[ngreyp-1].p[0] + 0.7 * lpnts[ngreyp-2].p[0]; + lpnts[ngreyp].v[0] = 0.3 * lpnts[ngreyp-1].v[0] + 0.7 * lpnts[ngreyp-2].v[0]; + lpnts[ngreyp++].w = 1.0; + + lpnts[ngreyp].p[0] = 0.7 * lpnts[ngreyp-2].p[0] + 0.3 * lpnts[ngreyp-3].p[0]; + lpnts[ngreyp].v[0] = 0.7 * lpnts[ngreyp-2].v[0] + 0.3 * lpnts[ngreyp-3].v[0]; + lpnts[ngreyp++].w = 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 cpll, cplv; /* Center point location and value */ - double kpwpos = 0.30; /* [0.30] White knee point location prop. towards center */ - double kpbpos = 0.15; /* [0.15] Black knee point location prop. towards center */ + 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 */ @@ -1069,20 +994,20 @@ glumknf = 1.0; printf("%sdbL = %f, dwL = %f\n", revrspl ? "(swapped) ": "", dbL,dwL); #endif - /* Center point location. Make lightly weighted */ - /* center the perceptual source, to try and maintain */ - /* the absolute source grey in the output, while */ - /* still allowing some of the knee compression/expansion to creep */ - /* into the other half. */ + /* 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; #ifdef PLOT_LMAP printf("cpll = %f, cplv = %f\n",cpll, cplv); #endif - /* Add weakish center point */ + +#ifdef NEVER /* Don't use a center point */ 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); @@ -1096,7 +1021,7 @@ glumknf = 1.0; // ~~ 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); - + //kwx = 0.0; //kbx = 0.0; //glumknf = 0.0; @@ -1110,6 +1035,7 @@ glumknf = 1.0; if (kbv < dbL) /* Sanity check */ kbv = dbL; + #ifdef PLOT_LMAP printf("using kbl = %f, kbv = %f\n",kbl, kbv); printf("using kwl = %f, kwv = %f\n",kwl, kwv); @@ -1119,110 +1045,68 @@ glumknf = 1.0; 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 */ + } - /* Create RSPL */ - { - datai il, ih; - datao ol, oh; - double avgdev[MXDO]; - int gres = 256; - - /* 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 */ - error("gamut: grey new_rspl failed"); - - il[0] = -1.0; /* Set possible input range */ - ih[0] = 101.0; - ol[0] = 0.0; /* Set normalisation output range */ - oh[0] = 100.0; + /* We now create the 1D rspl L map, that compresses or expands the luminence */ + /* range, independent of grey axis alignment, or gamut compression. */ + /* Because the rspl isn't symetrical when we swap X & Y, and we would */ + /* like a conversion from profile A to B to be the inverse of profile B to A */ + /* (as much as possible), we contrive here to always create a compression */ + /* RSPL, and create an inverse for it, and swap the two of them so that */ + /* the transform is correct and has an accurate inverse available. */ + { + datai il, ih; + datao ol, oh; + double avgdev[MXDO]; + int gres = 256; + + /* 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 */ + error("gamut: grey new_rspl failed"); + + il[0] = -1.0; /* Set possible input range */ + ih[0] = 101.0; + ol[0] = 0.0; /* Set normalisation output range */ + oh[0] = 100.0; #ifdef NEVER /* Dump out the L mapping points */ - { - int i; - printf("1D rspl L mapping points:\n"); - for (i = 0; i < ngreyp; i++) - printf("%d %f -> %f (w %f)\n",i,lpnts[i].p[0],lpnts[i].v[0],lpnts[i].w); - } + { + int i; + printf("1D rspl L mapping points:\n"); + for (i = 0; i < ngreyp; i++) + printf("%d %f -> %f (w %f)\n",i,lpnts[i].p[0],lpnts[i].v[0],lpnts[i].w); + } #endif - /* Create spline from the data points, with appropriate smoothness. */ - avgdev[0] = GAMMAP_RSPLAVGDEV; - if (s->grey->fit_rspl_w(s->grey, GAMMAP_RSPLFLAGS, lpnts, ngreyp, il, ih, &gres, ol, oh, 5.0, avgdev, NULL)) { - fprintf(stderr,"Warning: Grey axis mapping is non-monotonic - may not be very smooth ?\n"); - } - - /* Fine tune the rspl, to make sure that the white and black */ - /* point mapping is precise */ - { - co cp; - adjust1wb cx; /* Adjustment context */ - - /* Lookup actual black & white */ - cp.p[0] = sbL; - s->grey->interp(s->grey, &cp); - cx.awb[0] = cp.v[0]; - - cp.p[0] = swL; - s->grey->interp(s->grey, &cp); - cx.awb[1] = cp.v[0]; - - /* Set target black and white */ - cx.twb[0] = dbL; - cx.twb[1] = dwL; - - /* Fine tune the 3D->3D mapping */ - s->grey->re_set_rspl( - s->grey, /* this */ - 0, /* Combination of flags */ - (void *)&cx, /* Opaque function context */ - adjust1_wb_func /* Function to set from */ - ); - -#ifdef VERBOSE - if (verb) { - printf("Before tuning, L map White/Black is %f %f, should be %f %f\n", - cx.awb[1], cx.awb[0], dwL, dbL); - - /* Lookup fine tuned black & white */ - cp.p[0] = sbL; - s->grey->interp(s->grey, &cp); - cx.awb[0] = cp.v[0]; - - cp.p[0] = swL; - s->grey->interp(s->grey, &cp); - cx.awb[1] = cp.v[0]; - - printf("After tuning, L map White/Black is %f %f, should be %f %f\n", - cx.awb[1], cx.awb[0], dwL, dbL); - } -#endif /* VERBOSE */ - } - + /* Create spline from the data points, with appropriate smoothness. */ + avgdev[0] = GAMMAP_RSPLAVGDEV; + if (s->grey->fit_rspl_w(s->grey, GAMMAP_RSPLFLAGS, lpnts, ngreyp, il, ih, &gres, ol, oh, 5.0, avgdev, NULL)) { + fprintf(stderr,"Warning: Grey axis mapping is non-monotonic - may not be very smooth ?\n"); + } - /* Create an inverse mapping too, for reverse gamut and/or expansion. */ - il[0] = -1.0; /* Set possible input range */ - ih[0] = 101.0; - ol[0] = 0.0; /* Set normalisation output range */ - oh[0] = 100.0; + /* Create an inverse mapping too, for reverse gamut and/or expansion. */ + il[0] = -1.0; /* Set possible input range */ + ih[0] = 101.0; + ol[0] = 0.0; /* Set normalisation output range */ + oh[0] = 100.0; - if ((s->igrey = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) /* Allocate 1D -> 1D */ - error("gamut: igrey new_rspl failed"); - - /* 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 ((s->igrey = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) /* Allocate 1D -> 1D */ + error("gamut: igrey new_rspl failed"); + + /* 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 swapping of white and black points */ - rspl *tt = s->grey; - s->grey = s->igrey; - s->igrey = tt; - } + if (revrspl) { /* Swap to compensate for swapping of white and black points */ + rspl *tt = s->grey; + s->grey = s->igrey; + s->igrey = tt; } } @@ -1263,8 +1147,6 @@ glumknf = 1.0; if ((scl_gam = parttransgamut(s, sc_gam)) == NULL) { fprintf(stderr,"gamut map: parttransgamut failed\n"); - if (si_gam != sc_gam) - si_gam->del(si_gam); free(s); return NULL; } @@ -1275,8 +1157,6 @@ glumknf = 1.0; else { if ((sil_gam = parttransgamut(s, si_gam)) == NULL) { fprintf(stderr,"gamut map: parttransgamut failed\n"); - if (si_gam != sc_gam) - si_gam->del(si_gam); free(s); return NULL; } @@ -1289,6 +1169,7 @@ glumknf = 1.0; 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; datai il, ih; datao ol, oh; @@ -1297,8 +1178,7 @@ glumknf = 1.0; nearsmth *nsm = NULL; /* Returned list of near smooth points */ int nnsm; /* Number of near smoothed points */ double brad = 0.0; /* Black bend radius */ - gammapweights xpweights[14], xlpweights[14], xsweights[14]; - /* Explicit perceptial, lightnes pp. and sat. weights */ + gammapweights xpweights[14], xsweights[14]; /* Explicit perceptial and sat. weights */ gammapweights xwh[14]; /* Structure holding blended weights */ double smooth = 1.0; /* Level of 3D RSPL smoothing, blend of psmooth and ssmooth */ vrml *wrl = NULL; /* Gamut mapping illustration (hulls + guide vectors) */ @@ -1313,9 +1193,17 @@ typedef struct { #endif /* PLOT_3DKNEES */ /* Get the maximum number of points that will be created */ - nspts = near_smooth_np(NULL, scl_gam, sil_gam, d_gam, xvra, 4, surfpnts ? mapres : 0); + nspts = near_smooth_np(scl_gam, sil_gam, d_gam, xvra); - max_gpnts = nres + nspts; + rgridpts = 0; +#ifdef USE_BOUND + if (defrgrid >= 2) { + rgridpts = defrgrid * defrgrid * defrgrid + - (defrgrid -2) * (defrgrid -2) * (defrgrid -2); + } +#endif + + 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); @@ -1323,8 +1211,6 @@ typedef struct { if (sil_gam != scl_gam) sil_gam->del(sil_gam); scl_gam->del(scl_gam); - if (si_gam != sc_gam) - si_gam->del(si_gam); free(s); return NULL; } @@ -1365,10 +1251,8 @@ typedef struct { /* Create source grey axis point */ t = i/(nres - 1.0); -#ifdef NEVER /* Cover L = 0.0 to 100.0 */ t = ((100.0 * t) - sl_cs_bp[0])/(sl_cs_wp[0] - sl_cs_bp[0]); -#endif for (j = 0; j < 3; j++) gpnts[ngamp].p[j] = sl_cs_bp[j] + t * (sl_cs_wp[j] - sl_cs_bp[j]); @@ -1413,7 +1297,7 @@ typedef struct { gpnts[ngamp].w = wt; //printf("~1 t = %f, blended %f %f %f\n",t, gpnts[ngamp].v[0], gpnts[ngamp].v[1],gpnts[ngamp].v[2]); -#ifdef DUMP_GREY_AXIS_POINTS +#ifdef NEVER printf("Grey axis %d maps %f %f %f -> %f %f %f wit %f\n",ngamp, gpnts[ngamp].p[0], gpnts[ngamp].p[1], gpnts[ngamp].p[2], gpnts[ngamp].v[0], gpnts[ngamp].v[1], gpnts[ngamp].v[2], @@ -1486,7 +1370,6 @@ typedef struct { /* Convert from compact to explicit hextant weightings */ if (expand_weights(xpweights, pweights) - || expand_weights(xlpweights, lpweights) || expand_weights(xsweights, sweights)) { fprintf(stderr,"gamut map: expand_weights() failed\n"); s->grey->del(s->grey); @@ -1494,17 +1377,13 @@ typedef struct { if (sil_gam != scl_gam) sil_gam->del(sil_gam); scl_gam->del(scl_gam); - if (si_gam != sc_gam) - si_gam->del(si_gam); free(s); return NULL; } - - /* Create weights as blend between perceptual, lightness pp. and saturation */ - near_xwblend3(xwh, xpweights, gmi->gampwf, xlpweights, gmi->gamlpwf, - xsweights, gmi->gamswf); - if ((gmi->gampwf + gmi->gamlpwf + gmi->gamswf) > 0.1) - smooth = (gmi->gampwf * psmooth) + (gmi->gamlpwf * lpsmooth) + (gmi->gamswf * ssmooth); + /* Create weights as blend between perceptual and saturation */ + near_xwblend(xwh, xpweights, gmi->gampwf, xsweights, gmi->gamswf); + if ((gmi->gampwf + gmi->gamswf) > 0.1) + smooth = (gmi->gampwf * psmooth) + (gmi->gamswf * ssmooth); /* Tweak gamut mappings according to extra cmy cusp flags or rel override */ if (dst_cmymap != 0 || rel_oride != 0) { @@ -1516,7 +1395,7 @@ typedef struct { nsm = near_smooth(verb, &nnsm, scl_gam, sil_gam, d_gam, src_kbp, dst_kbp, dr_cs_bp, xwh, gmi->gamcknf, gmi->gamxknf, gmi->gamcpf > 1e-6, gmi->gamexf > 1e-6, - xvra, mapres, smooth, 1.10, surfpnts, il, ih, ol, oh); + xvra, mapres, smooth, il, ih, ol, oh); if (nsm == NULL) { fprintf(stderr,"Creating smoothed near points failed\n"); s->grey->del(s->grey); @@ -1524,11 +1403,58 @@ typedef struct { if (sil_gam != scl_gam) sil_gam->del(sil_gam); scl_gam->del(scl_gam); - if (si_gam != sc_gam) - si_gam->del(si_gam); free(s); return NULL; } + /* --------------------------- */ + + /* Make sure the input range to encompasss the guide vectors. */ + for (i = 0; i < nnsm; i++) { + for (j = 0; j < 3; j++) { + if (nsm[i].sv[j] < il[j]) + il[j] = nsm[i].sv[j];; + if (nsm[i].sv[j] > ih[j]) + ih[j] = nsm[i].sv[j]; + } + } + +#ifdef NEVER + if (verb) { + fprintf(stderr,"Input bounding box:\n"); + fprintfstderr,("%f -> %f, %f -> %f, %f -> %f\n", + il[0], ih[0], il[1], ih[1], il[2], ih[2]); + } +#endif + + /* 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. */ + { + 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; +#ifdef NEVER + if (verb) { + fprintf(stderr,"After incresing mapres to %d, input bounding box for 3D gamut mapping is:\n",mapres); + fprintf(stderr,"%f -> %f, %f -> %f, %f -> %f\n", + il[0], ih[0], il[1], ih[1], il[2], ih[2]); + } +#endif + } /* ---------------------------------------------------- */ /* Setup for diagnostic plot, that will have elements added */ @@ -1587,25 +1513,18 @@ typedef struct { for (i = 0; i < nnsm; i++) { double cpexf; /* The effective compression or expansion factor */ - /* Grid surface point */ - if (nsm[i].uflag != 0) { - cpexf = gmi->gamcpf; /* Assume compression */ + if (nsm[i].vflag == 0) { /* Unclear whether compression or expansion */ + /* Use larger to the the two factors */ + cpexf = gmi->gamcpf > gmi->gamexf ? gmi->gamcpf : gmi->gamexf; + + } else if (nsm[i].vflag == 1) { /* Compression */ + cpexf = gmi->gamcpf; + + } else if (nsm[i].vflag == 2) { /* Expansion */ + cpexf = gmi->gamexf; - /* Guide vector */ } else { - if (nsm[i].vflag == 0) { /* Unclear whether compression or expansion */ - /* Use larger to the the two factors */ - cpexf = gmi->gamcpf > gmi->gamexf ? gmi->gamcpf : gmi->gamexf; - - } else if (nsm[i].vflag == 1) { /* Compression */ - cpexf = gmi->gamcpf; - - } else if (nsm[i].vflag == 2) { /* Expansion */ - cpexf = gmi->gamexf; - - } else { - error("gammap: internal, unknown guide point flag"); - } + error("gammap: internal, unknown guide point flag"); } /* Compute destination value which is a blend */ @@ -1623,14 +1542,14 @@ typedef struct { gpnts[ngamp].p[j] = nsm[i].sv[j]; gpnts[ngamp].v[j] = nsm[i].div[j]; } - gpnts[ngamp++].w = nsm[i].w1; /* 1.01 for guide vectors, less for grid surface */ - + 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].uflag == 0 && nsm[i].vflag != 0) { /* Sub surface point is available */ + if (nsm[i].vflag != 0) { /* Sub surface point is available */ /* Compute destination value which is a blend */ /* between the source value and the knee adjusted destination */ @@ -1670,6 +1589,82 @@ typedef struct { 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 */ +#ifdef USE_BOUND + for (j = 0; j < 3; j++) { /* Set resolution for all axes */ + gres[j] = (mapres+1)/2; + avgdev[j] = GAMMAP_RSPLAVGDEV; + } + s->map = new_rspl(RSPL_NOFLAGS, 3, 3); /* Allocate 3D -> 3D */ + s->map->fit_rspl_w(s->map, GAMMAP_RSPLFLAGS, gpnts, ngamp, il, ih, gres, ol, oh, smooth, avgdev, NULL); + + /* Add input range grid surface anchor points to improve clipping behaviour. */ + if (defrgrid >= 2) { + DCOUNT(gc, 3, 3, 0, 0, defrgrid); + double cent[3]; + + sc_gam->getcent(d_gam, cent); + + DC_INIT(gc); + for (;;) { + /* If point is on the grid surface */ + if ( gc[0] == 0 || gc[0] == (defrgrid-1) + || gc[1] == 0 || gc[1] == (defrgrid-1) + || gc[2] == 0 || gc[2] == (defrgrid-1)) { + double grid2gamut, gamut2cent, ww; + co cp; + + /* Clip the point to the closest location on the source */ + /* colorspace gamut. */ + for (j = 0; j < 3; j++) + gpnts[ngamp].p[j] = il[j] + gc[j]/(defrgrid-1.0) * (ih[j] - il[j]); + sc_gam->nearest(sc_gam, cp.p, gpnts[ngamp].p); + + /* Then lookup the equivalent gamut mapped value */ + s->map->interp(s->map, &cp); + + for (j = 0; j < 3; j++) + gpnts[ngamp].v[j] = cp.v[j]; + + /* Compute the distance of the grid surface point to the to the */ + /* source colorspace gamut, as well as the distance from there */ + /* to the gamut center point. */ + for (grid2gamut = gamut2cent = 0.0, j = 0; j < 3; j++) { + double tt; + tt = gpnts[ngamp].p[j] - cp.p[j]; + grid2gamut += tt * tt; + tt = cp.p[j] - cent[j]; + gamut2cent += tt * tt; + } + grid2gamut = sqrt(grid2gamut); + gamut2cent = sqrt(gamut2cent); + + /* Make the weighting inversely related to distance, */ + /* to reduce influence on in gamut mapping shape, */ + /* while retaining some influence at the edge of the */ + /* grid. */ + ww = grid2gamut / gamut2cent; + if (ww > 1.0) + ww = 1.0; + + /* 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)) + break; + } + } +#else /* !USE_BOUND */ + printf("!!!! Warning - gammap boundary points disabled !!!!\n"); +#endif /* !USE_BOUND */ + /* --------------------------- */ /* Compute the output bounding values, and check input range hasn't changed */ for (i = 0; i < ngamp; i++) { @@ -1696,9 +1691,14 @@ typedef struct { #endif /* Create the final gamut mapping rspl. */ - /* [ How about converting to a delta filer ? ie. */ - /* create current filter, then create point list of delta from */ - /* smoothed value, filtering that and then un-deltering it ?? ] */ + /* [ The smoothing is not as useful as it should be, because */ + /* if it is increased it tends to push colors out of gamut */ + /* where they get clipped. Some cleverer scheme which makes */ + /* sure that smoothness errs on the side of more compression */ + /* is needed. - Addressed in nearsmth now ? ] */ + /* How about converting to a delta filer ? ie. */ + /* create curren filter, then create point list of delta from */ + /* smoothed value, filtering that and then un-deltering it ?? */ if (s->map != NULL) s->map->del(s->map); if (verb) @@ -1720,24 +1720,20 @@ typedef struct { /* (This isn't a good indication now that vectors have been adjusted */ /* to counteract the rspl smoothing at the edges.) */ if (verb) { - double de, avgde = 0.0, maxde = 0.0, num = 0.0; /* DE stats */ + double de, avgde = 0.0, maxde = 0.0; /* DE stats */ for (i = 0; i < nnsm; i++) { double av[3]; - if (nsm[i].uflag != 0) /* Ignore grid boundary points */ - continue; - /* Compute the mapping error */ dopartialmap2(s, av, nsm[i].sv); /* Just the rspl */ de = icmLabDE(nsm[i].div, av); avgde += de; - num++; if (de > maxde) maxde = de; } - printf("Gamut hull fit to guides: = avg %f, max %f\n",avgde/num,maxde); + printf("Gamut hull fit to guides: = avg %f, max %f\n",avgde/nnsm,maxde); } #endif /* CHECK_NEARMAP */ @@ -1834,10 +1830,6 @@ typedef struct { /* Show all neighbours */ wrl->start_line_set(wrl, 0); for (i = 0; i < nnsm; i++) { - - if (nsm[i].uflag != 0) /* Ignore grid boundary points */ - continue; - for (j = 0; j < XNNB; j++) { nearsmth *np = nsm[i].n[j]; /* Pointer to neighbor */ @@ -1873,8 +1865,6 @@ typedef struct { /* Locate the nearest source point */ for (ix = 0; ix < nnsm; ix++) { double dist = icmNorm33(pp, nsm[ix].sv); - if (nsm[i].uflag != 0) /* Ignore grid boundary points */ - continue; if (dist < bdist) { bdist = dist; bix = ix; @@ -1957,9 +1947,6 @@ typedef struct { double *ccc; double mdst[3]; - if (nsm[i].uflag != 0) /* Ignore grid boundary points */ - continue; - #if defined(SHOW_ACTUAL_VECTORS) || defined(SHOW_ACTUAL_VEC_DIFF) # ifdef SHOW_ACTUAL_VECTORS wrl->add_col_vertex(wrl, 0, nsm[i].sv, yellow); @@ -1996,8 +1983,6 @@ typedef struct { #if defined(SHOW_VECTOR_INDEXES) || defined(SHOW_SUB_PNTS) for (i = 0; i < nnsm; i++) { - if (nsm[i].uflag != 0) /* Ignore grid boundary points */ - continue; #ifdef SHOW_VECTOR_INDEXES { double cream[3] = { 0.7, 0.7, 0.5 }; @@ -2139,7 +2124,7 @@ typedef struct { icmMul3By3x4(vec, mat, vec); /* Intersect it with the source gamut */ - if (sil_gam->vector_isect(sil_gam, vec, cpoint, isect, + if (si_gam->vector_isect(si_gam, vec, cpoint, isect, NULL, NULL, NULL, NULL, NULL) == 0) { continue; } @@ -2194,7 +2179,7 @@ typedef struct { icmMul3By3x4(vec, mat, vec); /* Intersect it with the source gamut */ - if (sil_gam->vector_isect(sil_gam, vec, cpoint, isect, + if (si_gam->vector_isect(si_gam, vec, cpoint, isect, NULL, NULL, NULL, NULL, NULL) == 0) { warning("Ring %d vect %d diagnostic vector intersect failed",i,j); continue; @@ -2276,8 +2261,6 @@ typedef struct { if (sil_gam != scl_gam) sil_gam->del(sil_gam); scl_gam->del(scl_gam); - if (si_gam != sc_gam) - si_gam->del(si_gam); return s; } @@ -2468,20 +2451,6 @@ double *in /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Function to pass to rspl to re-set output values, */ -/* to adjust the 1D white and black points */ -static void -adjust1_wb_func( - void *pp, /* adjust1wb structure */ - double *out, /* output value to be adjusted */ - double *in /* corresponding input value */ -) { - adjust1wb *p = (adjust1wb *)pp; - - /* Do a linear re-mapping from actual to target */ - out[0] = (out[0] - p->awb[0]) * (p->twb[1] - p->twb[0])/(p->awb[1] - p->awb[0]) + p->twb[0]; -} - /* Function to pass to rspl to invert grey curve */ static void inv_grey_func( void *cntx, diff --git a/gamut/gamut.c b/gamut/gamut.c index 3fc1c97..dd4f43b 100644 --- a/gamut/gamut.c +++ b/gamut/gamut.c @@ -56,7 +56,7 @@ #define TRIANG_TOL 1e-10 /* [1e-10] Triangulation tollerance */ #define NORM_LOG_POW 0.25 /* [0.25] Normal, colorspace lopow value */ -#define RAST_LOG_POW 0.10 /* [0.10] Raster lopow value (is 0.05 too extreme ??) */ +#define RAST_LOG_POW 0.05 /* [0.05] Raster lopow value */ #undef TEST_CONVEX_HULL /* Use pure convex hull, not log hull */ @@ -156,7 +156,6 @@ static int getssvert(gamut *s, double *rad, double pos[3], double norm[3], int i static void startnexttri(gamut *s); static int getnexttri(gamut *s, int v[3]); static double volume(gamut *s); -static int write_to_vrml(gamut *s, vrml *wrl, double trans, int docusps); static int write_trans_vrml(gamut *s, char *filename, int doaxes, int docusps, void (*transform)(void *cntx, double out[3], double in[3]), void *cntx); static int write_vrml(gamut *s, char *filename, int doaxes, int docusps); @@ -175,9 +174,8 @@ static int compute_vector_isect(gamut *s, double *p1, double *p2, double *min, d static int compute_vector_isectns(gamut *s, double *p1, double *p2, gispnt *lp, int ll); static double log_scale(gamut *s, double ss); static int intersect(gamut *s, gamut *s1, gamut *s2); -static int nexpintersect(gamut *s, gamut *s1, gamut *s2); -static int expdstbysrcmdst(gamut *s, gamut *s1, gamut *s2, gamut *s3, - void (*cvect)(void *cntx, double *p2, double *p1), void *cntx); +static int compdstgamut(gamut *s, gamut *img, gamut *src, gamut *dst, int docomp, int doexp, + gamut *nedst, void (*cvect)(void *cntx, double *vec, double *pos), void *cntx); static int vect_intersect(gamut *s, double *rvp, double *ip, double *p1, double *p2, gtri *t); static void compgawb(gamut *s); @@ -658,8 +656,7 @@ int isRast /* Flag indicating Raster rather than colorspace */ s->getvert = getvert; s->volume = volume; s->intersect = intersect; - s->nexpintersect = nexpintersect; - s->expdstbysrcmdst = expdstbysrcmdst; + s->compdstgamut = compdstgamut; s->radial = radial; s->nradial = nradial; s->nearest = nearest; @@ -670,7 +667,6 @@ int isRast /* Flag indicating Raster rather than colorspace */ s->getwb = getwb; s->setcusps = setcusps; s->getcusps = getcusps; - s->write_to_vrml = write_to_vrml; s->write_vrml = write_vrml; s->write_trans_vrml = write_trans_vrml; s->write_gam = write_gam; @@ -1066,13 +1062,44 @@ double pp[3] /* rectangular coordinate of point */ } /* ------------------------------------ */ - -/* intersect implementation */ -/* Assumes s has been initialised. */ -static void intersect_imp(gamut *s, gamut *sa, gamut *sb) { +/* Initialise this gamut with the intersection of the */ +/* the two given gamuts. Return NZ on error. */ +/* Return 1 if gamuts are not compatible */ +/* (We assume that the this gamut is currently empty) */ +static int intersect(gamut *s, gamut *sa, gamut *sb) { int i, j, k; gamut *s1, *s2; + if (sa->compatible(sa, sb) == 0) + return 1; + + if IS_LIST_EMPTY(sa->tris) + triangulate(sa); + if IS_LIST_EMPTY(sb->tris) + triangulate(sb); + + s->isJab = sa->isJab; + + if (sa->isRast || sb->isRast) + s->isRast = 1; + + for (j = 0; j < 3; j++) + s->cent[j] = sa->cent[j]; + + /* Clear some flags */ + s->cswbset = 0; + s->cswbset = 0; + s->dcuspixs = 0; + + /* If either is a raster gamut, make it a raster gamut */ + if (s->isRast) + s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ + else + s->logpow = NORM_LOG_POW; + + /* Don't filter the points (gives a more accurate result ?) */ + s->nofilter = 1; + /* Add each source gamuts verticies that lie within */ /* the other gamut */ for (k = 0; k < 2; k++) { @@ -1092,7 +1119,6 @@ static void intersect_imp(gamut *s, gamut *sa, gamut *sb) { continue; pl = s2->nradial(s2, NULL, s1->verts[i]->p); - if (pl <= (1.0 + 1e-9)) { expand_gamut(s, s1->verts[i]->p); s1->verts[i]->f &= ~GVERT_ISOS; /* s1 vert is not outside s2 */ @@ -1135,142 +1161,6 @@ static void intersect_imp(gamut *s, gamut *sa, gamut *sb) { } END_FOR_ALL_ITEMS(tp1); } -} - -/* Initialise this gamut with the intersection of the */ -/* the two given gamuts. Return NZ on error. */ -/* Return 1 if gamuts are not compatible */ -/* (We assume that the this gamut is currently empty) */ -static int intersect(gamut *s, gamut *sa, gamut *sb) { - gamut *ss; - int j; - - if (sa->compatible(sa, sb) == 0) - return 1; - - if IS_LIST_EMPTY(sa->tris) - triangulate(sa); - if IS_LIST_EMPTY(sb->tris) - triangulate(sb); - - s->sres = sa->sres > sb->sres ? sa->sres : sb->sres; - - s->isJab = sa->isJab; - - /* Clear some flags */ - s->cswbset = 0; - s->cswbset = 0; - s->dcuspixs = 0; - - /* If either is a raster gamut, make it a raster gamut */ - if (sa->isRast || sb->isRast) - s->isRast = 1; - - if (s->isRast) { - s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ - s->no2pass = 1; /* Only do one pass */ - } else { - s->logpow = NORM_LOG_POW; /* Convex hull compression power */ - s->no2pass = 0; /* Do two passes */ - } - - for (j = 0; j < 3; j++) - s->cent[j] = sa->cent[j]; - - /* Grab white & black from whichever source has it */ - if (sb->cswbset) - ss = sb; - else if (sb->cswbset) - ss = sa; - - if (ss->cswbset) { - for (j = 0; j < 3; j++) { - s->cs_wp[j] = ss->cs_wp[j]; - s->cs_bp[j] = ss->cs_bp[j]; - s->cs_kp[j] = ss->cs_kp[j]; - } - s->cswbset = ss->cswbset; - } - - /* Don't filter the points (gives a more accurate result ?) */ - s->nofilter = 1; - - intersect_imp(s, sa, sb); - - if (sa->gawbset) { - compgawb(s); - } - - s->nofilter = 0; - - return 0; -} - -/* ------------------------------------ */ - -/* Initialise this gamut with neutral axis points from sa, */ -/* and then intersected with sb. */ -/* Return NZ on error. */ -/* Expand sb with neutral points, and then intersect */ -/* with sa. */ - -/* Return 1 if gamuts are not compatible */ -/* (We assume that the this gamut is currently empty) */ -static int nexpintersect(gamut *s, gamut *sa, gamut *sb) { - int i, j, k; - gamut *s1, *s2; - - if (sa->compatible(sa, sb) == 0) - return 1; - - if IS_LIST_EMPTY(sa->tris) - triangulate(sa); - if IS_LIST_EMPTY(sb->tris) - triangulate(sb); - - s->sres = sa->sres > sb->sres ? sa->sres : sb->sres; - - s->isJab = sa->isJab; - - /* If either is a raster gamut, make it a raster gamut */ - if (sa->isRast || sb->isRast) - s->isRast = 1; - - if (s->isRast) { - s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ - s->no2pass = 1; /* Only do one pass */ - } else { - s->logpow = NORM_LOG_POW; /* Convex hull compression power */ - s->no2pass = 0; /* Do two passes */ - } - - for (j = 0; j < 3; j++) - s->cent[j] = sa->cent[j]; - - /* Clear some flags */ - s->cswbset = 0; - s->cswbset = 0; - s->dcuspixs = 0; - - /* Don't filter the points (gives a more accurate result ?) */ - s->nofilter = 1; - - /* Number of points to generate */ - k = 10; - if (sa->nv <= 10) - k = 5; - - /* Generate points from black to white */ - for (i = 0; i < k; i++) { - double pp[3]; - double bf = i/(k-1.0); - - icmBlend3(pp, sa->cs_bp, sa->cs_wp, bf); - expand_gamut(s, pp); - } - - /* let intersect_imp to the hard work */ - intersect_imp(s, sa, sb); s->nofilter = 0; @@ -1279,36 +1169,49 @@ static int nexpintersect(gamut *s, gamut *sa, gamut *sb) { /* ------------------------------------ */ /* - Initialise this gamut with the image/destination gamut - expanded by the amount that dest colorspace is outside - the source colorspace gamut. + Initialise this gamut with a destination mapping gamut. + s1 is the image gamut (a possible sub-set of the src gamut) + s2 is the source gamut + s3 is the destination gamut + + if docomp: + this gamut will be set to the smaller of the img & dst gamuts + else + this gamut will be the img gamut + + if doexp + this gamut will be expanded by the amount dst is outside the src gamut. The vector direction of "inwards" is that returned by the - callback function as p1 -> p2 if it is supplied, or radially - inwards if it is not. p2 is a "center" point to compute depth to. + callback function if it is supplied, or radially inwards + if it is not. Return 1 if gamuts are not compatible. - (We assume that the _this_ gamut is currently empty) + (We assume that the this gamut is currently empty) */ #define MXNIS 40 /* Maximum raw intersections handled */ -static int expdstbysrcmdst( -gamut *s, /* Gamut to be expanded */ -gamut *s1, /* (Image) destination gamut to be expanded */ +static int compdstgamut( +gamut *s, /* This */ +gamut *s1, /* Image gamut, assumed to be a subset of source gamut */ gamut *s2, /* The source space gamut */ gamut *s3, /* The destination space gamut */ +int docomp, /* Flag, nz if compression is enabled */ +int doexp, /* Flag, nz if expansion is enabled. */ +gamut *nedst, /* Optional - if doexp, then expand nedst with non-expanded dst gamut */ void (*cvect)(void *cntx, double *p2, double *p1), /* Compression direction callback */ -void *cntx /* which returns p2 which is in desired direction from given p1 */ +void *cntx /* Returns p2 which is in desired direction from given p1 */ ) { +#ifdef STATS + int tedges, edgestested, testhits; +#endif int i, j, k; gamut *ss[3]; + gispnt lp1[MXNIS], lp2[MXNIS], lp3[MXNIS]; /* Intersection lists */ int ll1, ll2, ll3; /* Returned list length */ - gispnt lp1[MXNIS], lp2[MXNIS], lp3[MXNIS]; /* Lists of intersections */ int ii, jj, kk; /* List indexes */ -//printf("\n~1 expdstbysrcmdst() called\n"); - if (s1->compatible(s1, s2) == 0 || s1->compatible(s2, s3) == 0) return 1; @@ -1323,14 +1226,6 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */ s->isJab = s1->isJab; s->isRast = s1->isRast; - if (s->isRast) { - s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ - s->no2pass = 1; /* Only do one pass */ - } else { - s->logpow = NORM_LOG_POW; /* Convex hull compression power */ - s->no2pass = 0; /* Do two passes */ - } - for (j = 0; j < 3; j++) s->cent[j] = s1->cent[j]; @@ -1339,6 +1234,12 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */ s->cswbset = 0; s->dcuspixs = 0; + /* If s1 is a raster gamut, make output a raster gamut */ + if (s->isRast) + s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ + else + s->logpow = NORM_LOG_POW; + /* Don't filter the points (gives a more accurate result ?) */ s->nofilter = 1; @@ -1346,12 +1247,26 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */ ss[1] = s2; ss[2] = s3; - /* Use all the triangle vertices from the two/three gamuts */ +//printf("~1 compdstgamut docomp %d, doexp %d\n",docomp,doexp); + /* Reset the above/below flags */ + /* 1 = img above dst */ + /* 2 = img below dst */ + /* 4 = src above dst */ + /* 8 = src below dst */ + for (k = 0; k < 3; k++) { /* For img, src & dst */ + for (i = 0; i < ss[k]->nv; i++) { + if (!(ss[k]->verts[i]->f & GVERT_TRI)) + continue; + + ss[k]->verts[i]->as = 0; + } + } + + /* Use all the triangle vertices from the three gamuts */ /* as candidate points, because any of them might */ /* determine a surface feature. */ - for (k = 0; k < 3; k++) { + for (k = 0; k < 3; k++) { /* For img, src & dst */ - /* For each vertex */ for (i = 0; i < ss[k]->nv; i++) { double pp[3], ppv, p2[3]; double rr, r4; @@ -1361,39 +1276,22 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */ icmCpy3(pp, ss[k]->verts[i]->p); /* Point in question */ - if (k == 0) { /* Seed with all points from gamut */ - expand_gamut(s, pp); /* to ensure result can't be less. */ - } - //printf("\n~1 k %d, point %d: %f %f %f\n", k,i,pp[0],pp[1],pp[2]); - - /* Get the mapping vector */ if (cvect != NULL) - cvect(cntx, p2, pp); /* Get mapping direction to center */ + cvect(cntx, p2, pp); /* Get mapping direction */ else - icmCpy3(p2, ss[k]->cent); /* Radial vector to center */ - icmNormalize33(pp, pp, p2, 1.0); /* Make p2->pp length 1.0 */ - -//printf("~1 k %d, center %d: %f %f %f\n", k,i,p2[0],p2[1],p2[2]); - - /* Locate the intersecting segments for each gamut. */ - /* The returned parameter value will be >= 1.0 at and beyond the center point */ - /* and < 1.0 on the pp side. */ - /* We need intersections for all three for this to be a potential expansion point. */ - if ((ll1 = s1->vector_isectns(s1, pp, p2, lp1, MXNIS)) == 0) { /* Dest gamut */ -//printf("~1 no dst intersection\n"); - continue; - } - if ((ll2 = s2->vector_isectns(s2, pp, p2, lp2, MXNIS)) == 0) { /* Src space */ -//printf("~1 no sc intersection\n"); - continue; - } - - if ((ll3 = s3->vector_isectns(s3, pp, p2, lp3, MXNIS)) == 0) { /* Dst space */ -//printf("~1 no dc intersection\n"); - continue; - } - + icmCpy3(p2, ss[k]->cent); /* Radial vector */ + icmNormalize33(p2, p2, pp, 1.0); + + /* Locate the intersecting segments for each gamut */ + ll1 = ll2 = ll3 = 0; + ll1 = s1->vector_isectns(s1, pp, p2, lp1, MXNIS); /* Image gamut */ + if (doexp) /* For dst - src expansion */ + ll2 = s2->vector_isectns(s2, pp, p2, lp2, MXNIS); /* Src gamut */ + if (docomp || doexp) /* For img ^ dst or dst - src */ + ll3 = s3->vector_isectns(s3, pp, p2, lp3, MXNIS); /* Dst gamut */ + +//printf("~1 ll1 %d, ll2 %d, ll3 %d\n",ll1,ll2,ll3); #ifdef NEVER printf("img segments:\n"); for (ii = 0; ii < ll1; ii++) @@ -1405,120 +1303,264 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */ for (ii = 0; ii < ll3; ii++) printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp3[ii].pv,lp3[ii].dir,lp3[ii].edge,lp3[ii].tri->n); #endif + /* Now go through each image segment */ + for (ii = 0; ii < ll1; ii += 2) { - /* We're only interested in the most outside intersection of each surface */ - /* on the side of the point in question. (we're ignoring and complex */ - /* topology) */ - if (lp1[0].pv > (1.0 - 1e-8)) { -//printf("~1 dst point is on other side\n"); - continue; - } - if (lp2[0].pv > (1.0 - 1e-8)) { -//printf("~1 sc point is on other side\n"); - continue; - } - if (lp3[0].pv > (1.0 - 1e-8)) { -//printf("~1 sc point is on other side\n"); - continue; - } - - /* Make sure that sc is inside dc, */ - /* and sc is on or above dst */ - if (lp2[0].pv > (lp3[0].pv - 1e-8) - && lp2[0].pv <= (lp1[0].pv + 1e-8)) { - double ex[3], sf; - - /* Expand point by up to dst - src */ - icmSub3(ex, lp3[0].ip, lp2[0].ip); - - /* Make expansion proportional to how far dst */ - /* is to sc */ - sf = (1.0 - lp1[0].pv)/(1.0 - lp2[0].pv); - icmScale3(ex, ex, sf); - icmAdd3(pp, lp1[0].ip, ex); - -//printf("~1 expanding point by sf %f = %f %f %f\nb",sf,ex[0],ex[1],ex[2]); -//printf("~1 to %f %f %f\nb",pp[0],pp[1],pp[2]); + icmCpy3(pp, lp1[ii].ip); + ppv = lp1[ii].pv; + +//printf("~1 img pnt at %f\n",lp1[ii].pv); + if (docomp) { + + /* Locate a dst segment around or below img point */ + for (kk = 0; kk < ll3; kk += 2) { + if ((lp1[kk+1].pv + 1e-8) >= ppv) + break; + } + + if (kk >= ll3) { /* No dst segments or none below */ + ss[k]->verts[i]->as |= 1; /* img above dst */ +//printf("~1 img pnt is outside dst\n"); + continue; + } else { +//printf("~1 ing %f - %f, dst %f - %f\n", lp1[ii].pv,lp1[ii+1].pv,lp3[kk].pv,lp3[kk+1].pv); + /* Use the dst point if it is smaller */ + if (ppv < lp3[kk].pv) { + icmCpy3(pp, lp3[kk].ip); + ppv = lp3[kk].pv; +//printf("~1 dst is smaller %f\n",ppv); + ss[k]->verts[i]->as |= 1; /* img above dst */ + } else { +//printf("~1 ing is smaller %f\n",ppv); + ss[k]->verts[i]->as |= 2; /* img below dst */ + } + } + } + + if (nedst != NULL) + expand_gamut(nedst, pp); + + while (doexp) { + /* Locate a src segment that ing point lies in. */ + for (jj = 0; jj < ll2; jj += 2) { + if (lp1[ii].pv >= (lp2[jj].pv - 1e-8) + && lp1[ii].pv <= (lp2[jj+1].pv + 1e-8)) + break; + } + if (jj >= ll2) { + ss[k]->verts[i]->as |= 4; /* src above dst ?? */ + break; /* No overlapping segment */ + } + + /* Locate a dst segment that src point lies in */ + for (kk = 0; kk < ll3; kk += 2) { + if (lp2[jj].pv >= (lp3[kk].pv - 1e-8) + && lp2[jj].pv <= (lp3[kk+1].pv + 1e-8)) + break; + } + if (kk >= ll2) { + ss[k]->verts[i]->as |= 4; /* src above dst ?? */ + break; /* No overlapping segment */ + } + + /* if src is outside dst, do nothing */ + if (lp3[kk].pv >= lp2[jj].pv) { + ss[k]->verts[i]->as |= 4; /* src above dst */ + break; + } + + /* Expand point by dst - src */ + icmAdd3(pp, pp, lp3[kk].ip); + icmSub3(pp, pp, lp2[jj].ip); + ss[k]->verts[i]->as |= 8; /* src below dst */ +//printf("~1 expanding point by %f - %f\nb",lp3[kk].pv,lp2[jj].pv); + break; + } expand_gamut(s, pp); - } else { -//printf("~1 lp2 %f <= lp3 %f || lp2 %f > lp2 %f\n",lp2[0].pv,lp3[0].pv - 1e-8, lp2[0].pv, lp1[0].pv + 1e-8); -//printf("~1 dc is inside sc or sc is inside dst\n"); } } } - /* Generate points along the intersection of sc and dc, map them */ - /* to the image/dest gamut, and add them as well. This is to properly define */ - /* the edges of the expansion zone. */ +#ifdef STATS + tedges = edgestested = testhits = 0; +#endif - /* For sc on dc and then dc on sc */ - for (k = 0; k < 2; k++) { - gamut *ss1, *ss2; - gtri *tp1, *tp2; /* Triangle pointers */ + /* Add any intersection points between img/dst, and src/dst */ + for (k = 0; k < 4; k++) { + int mask; + gamut *sa, *sb; + gtri *tpa, *tpb; /* Triangle pointers */ if (k == 0) { - ss1 = s2; - ss2 = s3; - } else { - ss1 = s3; - ss2 = s2; + if (!docomp) + continue; + mask = 3; + sa = s1; /* img */ + sb = s3; /* dst */ + } else if (k == 1) { + if (!docomp) + continue; + mask = 3; + sa = s3; /* dst */ + sb = s1; /* img */ + } else if (k == 2) { + if (!doexp) + continue; + mask = 12; + sa = s2; /* src */ + sb = s3; /* dst */ + } else if (k == 3) { + if (!doexp) + continue; + mask = 12; + sa = s3; /* dst */ + sb = s2; /* src */ } - /* Find the edges that intersect the other gamut */ - tp1 = ss1->tris; - FOR_ALL_ITEMS(gtri, tp1) { /* For all ss1 triangles */ - - for (j = 0; j < 3; j++) { /* For all edges in ss1 triangle */ - /* If edge passes through the other gamut */ - if ((tp1->e[j]->v[0]->f ^ tp1->e[j]->v[1]->f) & GVERT_ISOS) { + /* Now find the edges that intersect the other gamut */ + tpa = sa->tris; + FOR_ALL_ITEMS(gtri, tpa) { - /* Exhaustive search of other triangles in ss2, */ - /* to find the one that the edge intersects with. */ - tp2 = ss2->tris; - FOR_ALL_ITEMS(gtri, tp2) { + for (j = 0; j < 3; j++) { /* For each edge */ +#ifdef STATS + tedges++; +#endif + /* If edge disposition flags aren't the same */ + if ((tpa->e[j]->v[0]->as ^ tpa->e[j]->v[1]->as) & mask + || (tpa->e[j]->v[0]->as & mask) == 3 || (tpa->e[j]->v[0]->as & mask) == 12 + || (tpa->e[j]->v[1]->as & mask) == 3 || (tpa->e[j]->v[1]->as & mask) == 12) { +#ifdef STATS + edgestested++; +#endif + /* Exhaustive search of other triangles */ + tpb = sb->tris; + FOR_ALL_ITEMS(gtri, tpb) { double pv; double pp[3]; /* Do a min/max intersection elimination test */ for (i = 0; i < 3; i++) { - if (tp2->mix[1][i] < tp1->mix[0][i] - || tp2->mix[0][i] > tp1->mix[1][i]) + if (tpb->mix[1][i] < tpa->mix[0][i] + || tpb->mix[0][i] > tpa->mix[1][i]) break; /* min/max don't overlap */ } if (i < 3) continue; /* Skip this triangle, it can't intersect */ - if (vect_intersect(ss1, &pv, pp, tp1->e[j]->v[0]->p, tp1->e[j]->v[1]->p, tp2) + if (vect_intersect(sa, &pv, pp, tpa->e[j]->v[0]->p, tpa->e[j]->v[1]->p, tpb) && pv >= (0.0 - 1e-10) && pv <= (1.0 + 1e-10)) { - double p2[3]; - - /* Got intersection point pp. */ - - /* Get the mapping vector */ + /* Process intersection point pp the same as the first loop */ + double ppv, p2[3]; + double rr, r4; +#ifdef STATS + testhits++; +#endif +//printf("\n~1 tri isxn point %f %f %f\n", pp[0],pp[1],pp[2]); if (cvect != NULL) - cvect(cntx, p2, pp); /* Get mapping direction to center */ + cvect(cntx, p2, pp); /* Get mapping direction */ else - icmCpy3(p2, ss[k]->cent); /* Radial vector to center */ - icmNormalize33(pp, pp, p2, 1.0); /* Make p2->pp length 1.0 */ - - /* Locate the intersecting segments for the img/dest gamut. */ - /* The returned parameter value will be >= 1.0 at and beyond */ - /* the center point and < 1.0 on the pp side. */ - /* We're only going to bother with the most outside point */ - if ((ll1 = s1->vector_isectns(s1, pp, p2, lp1, MXNIS)) == 0 - || lp1[0].pv > (1.0 - 1e-8)) { - continue; + icmCpy3(p2, sa->cent); /* Radial vector */ + icmNormalize33(p2, p2, pp, 1.0); + + /* Locate the intersecting segments for each gamut */ + ll1 = ll2 = ll3 = 0; + ll1 = s1->vector_isectns(s1, pp, p2, lp1, MXNIS); /* Image gamut */ + if (doexp) /* For dst - src expansion */ + ll2 = s2->vector_isectns(s2, pp, p2, lp2, MXNIS); /* Src gamut */ + if (docomp || doexp) /* For img ^ dst or dst - src */ + ll3 = s3->vector_isectns(s3, pp, p2, lp3, MXNIS); /* Dst gamut */ + +//printf("~1 ll1 %d, ll2 %d, ll3 %d\n",ll1,ll2,ll3); +#ifdef NEVER + printf("img segments:\n"); + for (ii = 0; ii < ll1; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp1[ii].pv,lp1[ii].dir,lp1[ii].edge,lp1[ii].tri->n); + printf("src segments:\n"); + for (ii = 0; ii < ll2; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp2[ii].pv,lp2[ii].dir,lp2[ii].edge,lp2[ii].tri->n); + printf("dst segments:\n"); + for (ii = 0; ii < ll3; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp3[ii].pv,lp3[ii].dir,lp3[ii].edge,lp3[ii].tri->n); +#endif + /* Now go through each image segment */ + for (ii = 0; ii < ll1; ii += 2) { + + icmCpy3(pp, lp1[ii].ip); + ppv = lp1[ii].pv; + +//printf("~1 img pnt at %f\n",lp1[ii].pv); + if (docomp) { + + /* Locate a dst segment around or below img point */ + for (kk = 0; kk < ll3; kk += 2) { + if ((lp1[kk+1].pv + 1e-8) >= ppv) + break; + } + + if (kk >= ll3) { /* No dst segments or none below */ +//printf("~1 img pnt is outside dst\n"); + continue; + } else { +//printf("~1 ing %f - %f, dst %f - %f\n", lp1[ii].pv,lp1[ii+1].pv,lp3[kk].pv,lp3[kk+1].pv); + /* Use the dst point if it is smaller */ + if (ppv < lp3[kk].pv) { + icmCpy3(pp, lp3[kk].ip); + ppv = lp3[kk].pv; +//printf("~1 dst is smaller %f\n",ppv); + } else { +//printf("~1 ing is smaller %f\n",ppv); + } + } + } + + if (nedst != NULL) + expand_gamut(nedst, pp); + + while (doexp) { + /* Locate a src segment that ing point lies in */ + for (jj = 0; jj < ll2; jj += 2) { + if (lp1[ii].pv >= (lp2[jj].pv - 1e-8) + && lp1[ii].pv <= (lp2[jj+1].pv + 1e-8)) + break; + } + if (jj >= ll2) { + break; /* No overlapping segment */ + } + + /* Locate a dst segment that src point lies in */ + for (kk = 0; kk < ll3; kk += 2) { + if (lp2[jj].pv >= (lp3[kk].pv - 1e-8) + && lp2[jj].pv <= (lp3[kk+1].pv + 1e-8)) + break; + } + if (kk >= ll2) { + break; /* No overlapping segment */ + } + + /* if src is outside dst, do nothing */ + if (lp3[kk].pv >= lp2[jj].pv) { + break; + } + + /* Expand point by dst - src */ + icmAdd3(pp, pp, lp3[kk].ip); + icmSub3(pp, pp, lp2[jj].ip); +//printf("~1 expanding point by %f - %f\nb",lp3[kk].pv,lp2[jj].pv); + break; + } + expand_gamut(s, pp); } - - expand_gamut(s, pp); } - } END_FOR_ALL_ITEMS(tp2); + } END_FOR_ALL_ITEMS(tpb); } } - } END_FOR_ALL_ITEMS(tp1); + } END_FOR_ALL_ITEMS(tpa); } +#ifdef STATS + printf("Total edges %d, edges tested %d, edge hits %d\n", tedges, edgestested, testhits); +#endif s->nofilter = 0; return 0; @@ -3907,7 +3949,7 @@ double *nin /* Normalised center relative point */ if (fabs(denom) < 1e-9) { /* Hmm. The ray is paralell to the triangle ? */ - error("radial_point: failed to intersect radial triangle, num %e, denom %e\n",num,denom); + error("radial_point: failed to intersect radial triangle\n"); } rv = num/denom; @@ -5582,7 +5624,7 @@ int ll /* Size of list. */ return 0; /* Too few to be useful */ } - /* Now we need to turn the raw intersections into sanitized segment pairs. */ + /* Now we need to turn th raw intersections into sanitized segment pairs. */ /* Sort the intersections by parameter value */ #define HEAP_COMPARE(A,B) (A.pv < B.pv) @@ -5759,67 +5801,7 @@ int ll /* Size of list. */ #endif /* INTERSECT_DEBUG */ /* ===================================================== */ - -/* Append gamut to an open VRML/X3d file */ -/* Return non-zero on error */ -static int write_to_vrml( - gamut *s, - vrml *wrl, - double trans, /* Transparency of gamut */ - int docusps /* Add cusps to vrml */ -) { - int i; - gtri *tp; /* Triangle pointer */ - - if IS_LIST_EMPTY(s->tris) - triangulate(s); - - if (docusps && s->cu_inited != 0) { - double ccolors[6][3] = { - { 1.0, 0.1, 0.1 }, /* Red */ - { 1.0, 1.0, 0.1 }, /* Yellow */ - { 0.1, 1.0, 0.1 }, /* Green */ - { 0.1, 1.0, 1.0 }, /* Cyan */ - { 0.1, 0.1, 1.0 }, /* Blue */ - { 1.0, 0.1, 1.0 } /* Magenta */ - }; - - for (i = 0; i < 6; i++) - wrl->add_marker(wrl, s->cusps[i], ccolors[i], 2.0); - } - - wrl->start_line_set(wrl, 0); - - /* Spit out the vertex values, in order. */ - for (i = 0; i < s->nv; i++) { - double out[3]; - - if (!(s->verts[i]->f & GVERT_TRI)) - continue; - - /* Show normal gamut surface */ - out[0] = s->verts[i]->p[0]; - out[1] = s->verts[i]->p[1]; - out[2] = s->verts[i]->p[2]; - wrl->add_vertex(wrl, 0, out); - } - - /* Add the surface triangles */ - tp = s->tris; - FOR_ALL_ITEMS(gtri, tp) { - int ix[3]; - ix[0] = tp->v[0]->tn; - ix[1] = tp->v[1]->tn; - ix[2] = tp->v[2]->tn; - wrl->add_triangle(wrl, 0, ix); - } END_FOR_ALL_ITEMS(tp); - - wrl->make_triangles_vc(wrl, 0, trans); - - return 0; -} - -/* Write gamut to a VRML/X3d file */ +/* Write to a VRML/X3d file */ /* Return non-zero on error */ static int write_vrml( gamut *s, @@ -5830,7 +5812,7 @@ int docusps /* Non-zero if cusp points are to be marked */ return write_trans_vrml(s, filename, doaxes, docusps, NULL, NULL); } -/* Write gamut to a VRML/X3d file */ +/* Write to a VRML/X3d file */ /* Return non-zero on error */ static int write_trans_vrml( gamut *s, @@ -5944,7 +5926,7 @@ void *cntx ix[2] = tp->v[2]->tn; wrl->add_triangle(wrl, 0, ix); } END_FOR_ALL_ITEMS(tp); -#endif /* !SHOW_BUCKETS */ +#endif /* SHOW_BUCKETS */ { double rgb[3]; diff --git a/gamut/gamut.h b/gamut/gamut.h index e70d898..467852f 100644 --- a/gamut/gamut.h +++ b/gamut/gamut.h @@ -36,8 +36,6 @@ #define MAXGAMN 10 /* Maximum gamut point neighbors returned */ #define NSLOTS 6 /* Number of maximum direction slots */ -struct _vrml *wrl; /* Declared in vrml.h, which may be #included after this */ - /* ------------------------------------ */ #define NODE_STRUCT \ int tag; /* Type of node, 1 = vertex, 2 = quad */ \ @@ -74,7 +72,7 @@ struct _gvert { double sp[3]; /* Point mapped to surface of unit sphere, relative to center */ double ch[3]; /* Point mapped for convex hull testing, relative to center */ - int as; /* Assert checking flag, expdstbysrcmdst flag */ + int as; /* Assert checking flag, compdstgamut flag */ }; typedef struct _gvert gvert; /* ------------------------------------ */ @@ -190,7 +188,7 @@ struct _gnn { struct _gispnt { double ip[3]; /* Intersecion Point */ double pv; /* Parameter value at intersection */ - int dir; /* Direction: 1 = into gamut, 0 = out of gamut */ + int dir; /* Direction: 1 = into gamut, 0 = out or gamut */ int edge; /* Edge: 2 = no isect, 1 = on edge, 0 = not on edge */ gtri *tri; /* Pointer to intersection triangle */ }; typedef struct _gispnt gispnt; @@ -315,14 +313,17 @@ struct _gamut { /* Initialise this gamut with the intersection of the */ /* the two given gamuts. */ - int (*nexpintersect)(struct _gamut *s, struct _gamut *s1, struct _gamut *s2); - /* Return s1 expanded with neutral axis points */ - /* and then intersected with s2. */ +#ifdef NEVER /* Deprecated */ + int (*expandbydiff)(struct _gamut *s, struct _gamut *s1, struct _gamut *s2, struct _gamut *s3, int docomp); + /* Initialise this gamut with a gamut which is s1 expanded */ + /* (but never reduced) by the distance from s2 to s3. */ + /* If docomp != 0, make gamut trace s3 if it's smaller than s1 */ +#endif - int (*expdstbysrcmdst)(struct _gamut *s, - struct _gamut *dst, struct _gamut *sc, struct _gamut *dc, - void (*cvect)(void *cntx, double *p2, double *p1), void *cntx); - /* Expand dst by ((dc - sc) > 0) */ + int (*compdstgamut)(struct _gamut *s, struct _gamut *img, struct _gamut *src, + struct _gamut *dst, int docomp, int doexpp, struct _gamut *nedst, + void (*cvect)(void *cntx, double *p2, double *p1), void *cntx); + /* Initialise this gamut with a destination mapping gamut. */ double (*radial)(struct _gamut *s, double out[3], double in[3]); /* return point on surface in same radial direction. */ @@ -383,8 +384,6 @@ struct _gamut { /* nz if no cusps available. */ /* Following return nz on error: */ - int (*write_to_vrml)(struct _gamut *s, struct _vrml *wrl, double trans, int docusps); - /* Append gamut surface to vrml. See also vrml->make_gamut_surface() etc. */ int (*write_vrml)(struct _gamut *s, char *filename, int doaxes, int docusps); /* Write to a VRML .wrl/.x3d file */ int (*write_gam)(struct _gamut *s, char *filename); /* Write to a CGATS .gam file */ diff --git a/gamut/maptest.c b/gamut/maptest.c index 433e719..3a33453 100644 --- a/gamut/maptest.c +++ b/gamut/maptest.c @@ -192,7 +192,6 @@ main(int argc, char *argv[]) { gmi.gampwf = 1.0; /* Gamut Perceptual Map weighting factor, 0.0 - 1.0 */ gmi.gamswf = 0.0; /* Gamut Saturation Map weighting factor, 0.0 - 1.0 */ } - gmi.gamlpwf = 0.0; /* Gamut Lightness preserving perceptual Map whtg. factor, 0.0 - 1.0 */ gmi.satenh = 0.0; /* Saturation enhancement factor */ gmi.desc = "mapetest"; gmi.icci = 0; @@ -206,7 +205,7 @@ main(int argc, char *argv[]) { 0, 0, /* Normal black points */ 0, /* Normal CMY cusp mapping */ 0, /* No relative weighting override */ - 19, /* rspl resolution of 17 */ + 17, /* rspl resolution of 17 */ NULL, /* No input range override */ NULL, gammapwrl /* Diagnostic plot */ diff --git a/gamut/nearsmth.c b/gamut/nearsmth.c index c65704e..c0bd2be 100644 --- a/gamut/nearsmth.c +++ b/gamut/nearsmth.c @@ -52,7 +52,6 @@ #include <fcntl.h> #include <string.h> #include <math.h> -#include "counters.h" #include "icc.h" #include "numlib.h" #include "rspl.h" @@ -61,14 +60,12 @@ #include "vrml.h" #undef SAVE_VRMLS /* [Und] Save various vrml's */ -#undef PLOT_SMOOTHING_CHANGE /* [Und] Dest point change due to smoothing in "dst_smvec.wrl" */ #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 /* [Und] Show the neighborhood point group in src */ #undef SHOW_NEIGB_WEIGHTS /* [Und] Show the weighting for each point of neighbours in turn */ #undef DIAG_POINTS /* [Und] Short circuite mapping and show vectors of various */ @@ -81,8 +78,9 @@ #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 VECADJPASSES 8 /* [8] Vector smoothing and adjust passes. */ -#define RSPLPASSES 4 /* [4] Number of rspl smoothing & adjustment passes */ +#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 smoothing to aim for */ #define SHRINK 5.0 /* [5.0] Shrunk destination evect surface factor */ #define CYLIN_SUBVEC /* [Def] Make sub-vectors always cylindrical direction */ @@ -116,7 +114,7 @@ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ #if defined(SAVE_VRMLS) && defined(PLOT_MAPPING_INFLUENCE) -static void create_influence_plot(nearsmth *smp, int nmpts, int mapres); +static void create_influence_plot(nearsmth *smp, int nmpts); #endif /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ @@ -196,7 +194,7 @@ double sumpow /* Sum power. 0.0 == 2.0 */ /* Compute the LCh differences squared of in1 - in2 */ /* (This is like the CIE DE94) */ -static void diffLChsq( +static void diffLCh( double out[3], double in1[3], /* Destination location */ double in2[3] /* Source location */ @@ -260,7 +258,7 @@ double dcratio, /* Depth compression ratio of mapping */ double dxratio /* Depth expansion ratio of mapping */ ) { double a_o; - double va, vr, vd, vv = 0.0; + 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, */ @@ -269,7 +267,7 @@ double dxratio /* Depth expansion ratio of mapping */ 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 */ - vr = wdesq(dtp, drv, w->rl.l, w->rl.c, w->rl.h, SUM_POW); + 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 @@ -278,16 +276,17 @@ double dxratio /* Depth expansion ratio of mapping */ /* Diagnostic values */ p->dbgv[0] = va; p->dbgv[1] = vr; - p->dbgv[2] = vd; + p->dbgv[2] = vl; + p->dbgv[3] = vd; - vv = va + vr + vd; /* Sum of squares */ -// vv = sqrt(va) + sqrt(vr) + sqrt(vd); /* Linear sum is better ? */ -// vv = pow(va, 0.7) + pow(vr, 0.7) + pow(vd, 0.7); /* Linear sum is better ? */ + 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 vr = %f from drv %f %f %f, weights %f %f %f\n", vr, drv[0], drv[1], drv[2], w->rl.l, w->rl.c, w->rl.h); + 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 */ @@ -301,13 +300,12 @@ double dxratio /* Depth expansion ratio of mapping */ /* and cusp mapping function. */ struct _smthopt { /* optimisation */ - int debug; /* debug flag */ int pass; /* Itteration round */ int ix; /* Index of point being optimized */ nearsmth *p; /* Point being optimised */ int useexp; /* Flag indicating whether expansion is permitted */ - double *wn; /* Target of weighted nearest */ - gamut *wngam; /* for optfunc1 and optfunc1a */ + int debug; /* debug flag */ + gamut *shgam; /* for optfunc1a */ /* Setup state */ int isJab; /* Flag indicating Jab rather than Lab space */ @@ -336,7 +334,7 @@ struct _smthopt { }; typedef struct _smthopt smthopt; -static void init_ce(smthopt *s, gamut *sc_gam, gamut *si_gam, gamut *d_gam, int src_kbp, int dst_kbp, double d_bp[3]); +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]); @@ -359,97 +357,84 @@ static void spow3(double *out, double *in, double ex) { } } -/* Absolute error function, used by optfunc1() & optfunc1a() */ -static double aerrf( - nearsmth *p, - double *dv, - double *sv +/* 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 ) { - double delch[3], rv; + 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 source gamut value */ + comp_ce(s, ddv, ddv, &p->wt); +//printf("~1 after cusp rot 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(dv, sv, p->wt.ra.l, p->wt.ra.c, p->wt.ra.h, SUM_POW); + 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; /* Extra power when L de is over thr */ + double ppp = p->wt.a.lxpow; double thr = p->wt.a.lxthr; /* Xover between normal and power */ double sumpow = SUM_POW; - double del; - diffLChsq(delch, dv, sv); - del = sqrt(delch[0]); /* delta L */ + 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], 1.0 + (ppp - 1.0) * del/(del + thr)); + 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], 1.0 + (ppp - 1.0) * del/(del + thr)) + 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], (1.0 + (ppp - 1.0) * del/(del + thr)) * sumpow) + + 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 - return rv; -} - -/* Powell optimisation function for setting minimal absolute error target point, */ -/* with a correction for swap. */ -/* We get a 2D plane in the 3D space, of the destination point, */ -/* who's location we are optimizing to wngam. */ -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 wngam surface */ - 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]); - s->wngam->radial(s->wngam, 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 source gamut value */ - comp_ce(s, ddv, ddv, &p->wt); -//printf("~1 after cusp rot got %f %f %f\n",ddv[0],ddv[1],ddv[2]); - } - - rv = aerrf(p, ddv, s->wn); if (s->debug) - printf("debug: rv = %f from %f %f %f -> %f %f %f\n",rv, s->wn[0], s->wn[1], s->wn[2], ddv[0], ddv[1], ddv[2]); + 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->wm[0], p->wm[1], p->wm[2], 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, */ -/* with no correction for swap. */ +/* 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 to wngam. */ +/* who's location we are optimizing. */ static double optfunc1a( void *fdata, double *_dv @@ -458,8 +443,9 @@ double *_dv 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 wngam surface */ - double rv; /* Out of gamut, return value */ + 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]; @@ -467,16 +453,48 @@ double *_dv dv[0] = 50.0; icmMul3By3x4(dv, p->m3d, dv); -//if (s->debug) printf("~1 optfunc1a got 2D %f %f -> 3D %f %f %f\n", _dv[0], _dv[1], dv[0], dv[1], dv[2]); - s->wngam->radial(s->wngam, ddv, dv); /* Map to shgam surface to check current location */ -//if (s->debug) printf("~1 optfunc1a got %f %f %f -> surface %f %f %f\n", dv[0], dv[1], dv[2], ddv[0], ddv[1], ddv[2]); +//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); - rv = aerrf(p, ddv, s->wn); + 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 -> %f %f %f\n",rv, s->wn[0], s->wn[1], s->wn[2], ddv[0], ddv[1], ddv[2]); + printf("debug: rv = %f from %f %f %f\n",rv, ddv[0], ddv[1], ddv[2]); -//if (s->debug) printf("~1 sv %4.2f %4.2f %4.2f, ddv %4.2f %4.2f %4.2f\n", p->wm[0], p->wm[1], p->wm[2], 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; } @@ -545,7 +563,7 @@ double *dv /* 3D Location being evaluated */ } } -/* Powell optimisation function for overall non-relative smoothed error optimization. */ +/* Powell optimisation function for non-relative error optimization. */ /* We get a 2D point in the 3D space. */ static double optfunc2( void *fdata, @@ -564,8 +582,8 @@ double *_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 */ @@ -579,13 +597,12 @@ double *_dv 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("~1 va = %f, vr = %f, vd = %f\n", p->dbgv[0], p->dbgv[1], p->dbgv[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); +//printf("~1 rv = %f\n\n",rv); return rv; } @@ -595,7 +612,6 @@ double *_dv static void init_ce( smthopt *s, /* Context for cusp mapping being set. */ gamut *sc_gam, /* Source colorspace gamut */ -gamut *si_gam, /* Source image 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 */ @@ -624,9 +640,9 @@ double d_bp[3] /* Override destination target black point (may be NULL) */ /* Set some default values for src white/black/grey */ - /* Get the colorspace white and black point info */ + /* Get the white and black point info */ if (src_kbp) { - if (sc_gam->getwb(sc_gam, s->cusps[0][6], NULL, s->cusps[0][7], NULL, NULL, NULL) != 0) { + 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; @@ -634,7 +650,7 @@ double d_bp[3] /* Override destination target black point (may be NULL) */ s->donaxis = 0; } } else { - if (sc_gam->getwb(sc_gam, s->cusps[0][6], s->cusps[0][7], NULL, NULL, NULL, NULL) != 0) { + 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; @@ -644,7 +660,7 @@ double d_bp[3] /* Override destination target black point (may be NULL) */ } if (dst_kbp) { - if (d_gam->getwb(d_gam, s->cusps[1][6], NULL, s->cusps[1][7], NULL, NULL, NULL) != 0) { + 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; @@ -652,7 +668,7 @@ double d_bp[3] /* Override destination target black point (may be NULL) */ s->donaxis = 0; } } else { - if (d_gam->getwb(d_gam, s->cusps[1][6], s->cusps[1][7], NULL, NULL, NULL, NULL) != 0) { + 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; @@ -664,21 +680,6 @@ double d_bp[3] /* Override destination target black point (may be NULL) */ icmCpy3(s->cusps[1][7], d_bp); } -#ifdef NEVER - { - double iwp[3] = { -1, -1, -1}, ibp[3] = { -1, -1, -1}; - if (src_kbp) { - si_gam->getwb(si_gam, NULL, NULL, NULL, iwp, NULL, ibp); - } else { - si_gam->getwb(si_gam, NULL, NULL, NULL, iwp, ibp, NULL); - } - printf("~1 src white = %f, black = %f\n",s->cusps[0][6][0],s->cusps[0][7][0]); - printf("~1 img white = %f, black = %f\n",s->cusps[0][6][0],s->cusps[0][7][0]); - printf("~1 dst white = %f, black = %f\n",s->cusps[1][6][0],s->cusps[1][7][0]); - } - -#endif /* NEVER */ - /* 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; @@ -1165,7 +1166,6 @@ gammapweights *src NSCOPY(r.rdl); NSCOPY(r.rdh); - NSCOPY(r.dsm); NSCOPY(d.co); NSCOPY(d.xo); @@ -1207,50 +1207,6 @@ gammapweights *src2, double wgt2 NSBLEND(r.rdl); NSBLEND(r.rdh); - NSBLEND(r.dsm); - - NSBLEND(d.co); - NSBLEND(d.xo); - - NSBLEND(f.x); -#undef NSBLEND -} - -/* Blend a three groups of individual weights into one, given three weightings */ -void near_wblend3( -gammapweights *dst, -gammapweights *src1, double wgt1, -gammapweights *src2, double wgt2, -gammapweights *src3, double wgt3 -) { - -#define NSBLEND(xxx) dst->xxx = wgt1 * src1->xxx + wgt2 * src2->xxx + wgt3 * src3->xxx - - NSBLEND(c.w.l); - NSBLEND(c.w.c); - NSBLEND(c.w.h); - NSBLEND(c.tw); - 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(r.dsm); NSBLEND(d.co); NSBLEND(d.xo); @@ -1393,9 +1349,8 @@ void tweak_weights(gammapweights out[14], int dst_cmymap, int rel_oride) { } if (rel_oride == 1) { /* A high saturation "clip" like mapping */ - out[i].r.rdl = 1.0; /* No relative neighbourhood/smoothing */ - out[i].r.rdh = 1.0; /* No relative neighbourhood/smoothing */ - out[i].r.dsm = 0.0; /* No relative neighbourhood/smoothing */ + 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 */ @@ -1406,7 +1361,7 @@ void tweak_weights(gammapweights out[14], int dst_cmymap, int rel_oride) { } } -/* Blend two expanded groups of individual weights into one */ +/* Blend a two expanded groups of individual weights into one */ void near_xwblend( gammapweights *dst, gammapweights *src1, double wgt1, @@ -1417,18 +1372,6 @@ gammapweights *src2, double wgt2 near_wblend(&dst[i], &src1[i], wgt1, &src2[i], wgt2); } -/* Blend three expanded groups of individual weights into one */ -void near_xwblend3( -gammapweights *dst, -gammapweights *src1, double wgt1, -gammapweights *src2, double wgt2, -gammapweights *src3, double wgt3 -) { - int i; - for (i = 0; i < 14; i++) - near_wblend3(&dst[i], &src1[i], wgt1, &src2[i], wgt2, &src3[i], wgt3); -} - /* Convert overall, hue dom & l dom to iweight */ static void comp_iweight(iweight *iw, double o, double h, double l) { double c, lc; @@ -1593,15 +1536,14 @@ void interp_xweights(gamut *gam, gammapweights *out, double pos[3], } } -/* Callback used by expdstbysrcmdst() to establish the expected compression */ -/* mapping direction. p2 should be the center point, so depth from the center */ -/* can be computed. We return a point on the neutral axis. */ +/* 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, lv[3]; + double vv, gv[3], lv[3]; smthopt *s = (smthopt *)cntx; gammapweights out; @@ -1617,26 +1559,18 @@ double *p1 /* Given point */ /* 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); - /* Normalise l * c weight to sum to 1.0 */ - vv = fabs(out.ra.l + out.ra.c); - if (vv < 1e-7) { /* Hmm. */ - out.ra.l = out.ra.c = 0.5; - } else { - out.ra.l /= vv; - out.ra.c /= vv; - } + icmSub3(gv, s->cusps[0][8], p1); /* Grey vector */ + icmNormalize3(gv, gv, out.ra.c); - /* Make p2 the weighted sum of equivalent L value and grey value on */ - /* the neutral axis. */ - icmScale3(lv, lv, out.ra.l); - icmScale3(p2, s->cusps[0][8], out.ra.c); - icmAdd3(p2, p2, lv); + icmAdd3(p2, gv, p1); + icmAdd3(p2, lv, p2); /* Combined destination */ //printf("~1 p2 %f %f %f\n", p2[0], p2[1], p2[2]); } @@ -1753,23 +1687,18 @@ double p1[3] /* Point */ /* ============================================ */ /* Return the maximum number of points that will be generated */ -/* (This isn't accurate due to manipulation of the gamuts in nearsmth!) */ int near_smooth_np( - gamut **pp_gam, /* Return gamut that was used for points */ gamut *sc_gam, /* Source colorspace gamut */ gamut *si_gam, /* Source image gamut (== sc_gam if none) */ - gamut *dc_gam, /* Destination colorspace gamut */ - double xvra, /* Extra vertex ratio */ - int gmult, /* Guide point multiplier, typically 4 */ - int surfgres /* surface grid point resolution, 0 for 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; - int hsurfgres = (surfgres + 1)/2; /* near_smooth uses half */ nspts = sc_gam->nverts(sc_gam); nipts = si_gam->nverts(si_gam); - ndpts = dc_gam->nverts(dc_gam); + ndpts = d_gam->nverts(d_gam); p_gam = sc_gam; /* Target number of points is max of any gamut */ @@ -1785,16 +1714,6 @@ int near_smooth_np( xvra = ntpts/(double)nspts; nmpts = p_gam->nssverts(p_gam, xvra); /* Stratified Sampling source points */ - nmpts *= gmult; /* Allow for sub-surface points etc. */ - - if (hsurfgres >= 4) { - nmpts += hsurfgres * hsurfgres * hsurfgres - - (hsurfgres -4) * (hsurfgres -4) * (hsurfgres -4); - } - - if (pp_gam != NULL) - *pp_gam = p_gam; - return nmpts; } @@ -1807,7 +1726,7 @@ 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 *dc_gam, /* Destination 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 */ @@ -1819,22 +1738,22 @@ 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 */ -double gexp, /* Grid expansion ratio, none = 1.0 */ -int surfpnts, /* Flag - add surface grid points */ -datai map_il, /* Return expanded input range */ +datai map_il, /* Preliminary rspl input range */ datai map_ih, -datao map_ol, /* Return expanded output range */ +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 *src_gam; /* Intersection of src and img gamut gamut */ - gamut *dst_gam; /* Modified destination gamut suitable for mapping from src_gam. */ - /* If compression, this is the intersection of src_gam and dc_gam. */ - /* If expansion, this is the src_gam expanded by dc_gam - sc_gam. */ - gamut *nedst_gam;/* Same as above, but not expanded. */ - int mxnmpts; /* Allocated number of mapping points */ + 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; @@ -1843,21 +1762,39 @@ datao map_oh double codf; /* Itteration overshoot/damping factor */ double mxmv; /* Maximum a point gets moved */ int nmxmv; /* Number of maxmoves less than stopping threshold */ - int dmapres = 1; /* Change in mapres when applying gexp */ - int hmapres; /* Half mapres */ - int hdmapres; /* Half change in mapres */ - rspl *lastmap = NULL; /* Last gamut mapping map created, if any */ /* Check gamuts are compatible */ - if (sc_gam->compatible(sc_gam, dc_gam) == 0 + 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; } - mxnmpts = near_smooth_np(&p_gam, sc_gam, si_gam, dc_gam, xvra, 1, surfpnts ? mapres : 0); - nmpts = 0; + { + 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 ? */ @@ -1868,106 +1805,80 @@ datao map_oh /* 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, si_gam, dc_gam, src_kbp, dst_kbp, d_bp); + init_ce(&opts, sc_gam, d_gam, src_kbp, dst_kbp, d_bp); /* Allocate our guide points */ - if ((smp = (nearsmth *)calloc(mxnmpts, sizeof(nearsmth))) == NULL) { + 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 image gamut intersected */ - /* with the source colorspace gamut, in case something strange with the */ - /* image gamut. (gammap.c may have already done this) */ - src_gam = sc_gam; /* Alias to source space gamut */ + /* 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 ((src_gam = new_gamut(0.0, 0, 0)) == NULL) { + if ((sci_gam = new_gamut(0.0, 0, 0)) == NULL) { fprintf(stderr,"gamut map: new_gamut failed\n"); free_nearsmth(smp, nmpts); *npp = 0; return NULL; } - src_gam->intersect(src_gam, si_gam, sc_gam); + sci_gam->intersect(sci_gam, sc_gam, si_gam); #ifdef SAVE_VRMLS { - char src_gam_name[40] = "si_gam"; - - printf("###### gamut/nearsmth.c: writing diagnostic si_gam%s, src_gam%s\n",vrml_ext(),vrml_ext()); - - strcat(src_gam_name, vrml_ext()); - src_gam->write_vrml(si_gam, src_gam_name, 1, 0); - - strcpy(src_gam_name, "src_gam"); - strcat(src_gam_name, vrml_ext()); - src_gam->write_vrml(src_gam, src_gam_name, 1, 0); + char sci_gam_name[40] = "sci_gam"; + strcat(sci_gam_name, vrml_ext()); + printf("###### gamut/nearsmth.c: writing diagnostic sci_gam%s and di_gam%s\n",vrml_ext(),vrml_ext()); + sci_gam->write_vrml(sci_gam, sci_gam_name, 1, 0); } #endif } - dst_gam = src_gam; /* Default no compress or expand */ - - /* non-expanded dst_gam for testing double back img points against: */ - nedst_gam = src_gam; /* Default same as dst_gam */ - - /* Convert dst_gam to compress and/or expand target for mapping src_gam to. */ + di_gam = sci_gam; /* Default no compress or expand */ if (usecomp || useexp) { - - if ((nedst_gam = dst_gam = new_gamut(0.0, 0, 0)) == NULL) { + if ((di_gam = new_gamut(0.0, 0, 0)) == NULL) { fprintf(stderr,"gamut map: new_gamut failed\n"); - if (src_gam != sc_gam) - src_gam->del(src_gam); + if (si_gam != sc_gam) + sci_gam->del(sci_gam); free_nearsmth(smp, nmpts); *npp = 0; return NULL; } - /* For compression only, nedst_gam and dst_gam are smaller of src_gam and dc_gam space. */ - /* Augment the dst_gam with neutral axis points in case the source gamut */ - /* has a "spike" that separates it from the neutral axis, allowing */ - /* mapping. */ - nedst_gam->nexpintersect(nedst_gam, dc_gam, src_gam); - if (useexp) { - /* No image gamut - dest colorspace is target */ - if (si_gam == sc_gam) { - dst_gam = dc_gam; /* Expanded dest is colorspace dest */ - - /* There is an image gamut, so */ - /* Expand nedst_gam to create dst_gam expanded in proportion to where */ - /* dc_gam is outside sc_gam */ - } else { - if ((dst_gam = new_gamut(0.0, 0, 0)) == NULL) { - fprintf(stderr,"gamut map: new_gamut failed\n"); - if (src_gam != sc_gam) - src_gam->del(src_gam); - free_nearsmth(smp, nmpts); - *npp = 0; - return NULL; - } - - /* Initialise this gamut with the nedst_gam expanded by ((dc_gam - sc_gam) > 0) */ - dst_gam->expdstbysrcmdst(dst_gam, nedst_gam, sc_gam, dc_gam, cvect, &opts); + /* 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); + free_nearsmth(smp, nmpts); + *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 { - char dst_gam_name[30] = "dst_gam"; - - printf("###### gamut/nearsmth.c: writing diagnostic dst_gam%s, nedst_gam%s\n",vrml_ext(),vrml_ext()); - strcat(dst_gam_name, vrml_ext()); - dst_gam->write_vrml(dst_gam, dst_gam_name, 1, 0); - - strcpy(dst_gam_name, "nedst_gam"); - strcat(dst_gam_name, vrml_ext()); - nedst_gam->write_vrml(nedst_gam, dst_gam_name, 1, 0); + char di_gam_name[30] = "di_gam"; + strcat(di_gam_name, vrml_ext()); + di_gam->write_vrml(di_gam, di_gam_name, 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 < mxnmpts; i++) { + 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 */ @@ -1979,12 +1890,12 @@ datao map_oh //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 = src_gam->radial(src_gam, imv, imv); + imr = sci_gam->radial(sci_gam, imv, imv); } /* If point is within non-expanded modified destination gamut, */ - /* then it is a "double back"/convex image point, and should be ignored. */ - if (nedst_gam->radial(nedst_gam, NULL, imv) > (imr + 1e-4)) { + /* 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; @@ -1992,18 +1903,17 @@ datao map_oh /* Lookup radialy equivalent point on modified destination gamut, */ /* in case we need it for compression or expansion */ - smp[i].drr = dst_gam->radial(dst_gam, smp[i].drv, imv); + 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].uflag = smp[i].vflag = smp[i].gflag = 0; + 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].w1 = 1.0; - smp[i].sgam = src_gam; - smp[i].dgam = src_gam; - smp[i].dcgam = dc_gam; + smp[i].sgam = sci_gam; + smp[i].dgam = sci_gam; + smp[i].mapres = mapres; VB(("In Src %d = %f %f %f\n",i,smp[i].sv[0],smp[i].sv[1],smp[i].sv[2])); @@ -2017,7 +1927,7 @@ datao map_oh double mv[3], ml; /* Radial inward mapping vector */ double dir; - icmSub3(mv, src_gam->cent, smp[i].sv); /* Vector to center */ + icmSub3(mv, sci_gam->cent, smp[i].sv); /* Vector to center */ ml = icmNorm3(mv); /* It's length */ if (ml > 0.001) { @@ -2037,7 +1947,6 @@ datao map_oh 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]; - smp[i].w1 = 1.01; /* Use 1.01 as marker value */ 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])); @@ -2046,17 +1955,17 @@ datao map_oh *npp = nmpts; /* Don't need this anymore */ - if (nedst_gam != src_gam && nedst_gam != dst_gam) - nedst_gam->del(nedst_gam); - nedst_gam = NULL; + 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 (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); + 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; } @@ -2069,15 +1978,6 @@ datao map_oh interp_xweights(opts.sgam, &smp[i].wt, smp[i]._sv, opts.xwh, &opts, 0); } - /* ~~ would be nice to eliminate the need for dst_gam that is the intersection - * of dc_gam and sc/img_gam here. Problem is determining expansion vector - * direction in a way that is consistent with the absolute error weighting. - * - * For the moment leave the current appoach of using the dst_gam that has been - * expanded in proportion to dc_gam - sc_gam in cvec() direction, since - * the absolute error weighting is use to map the sv to that surface. - */ - 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")); @@ -2100,7 +2000,7 @@ datao map_oh /* Compute the cusp rotated version of the cspace/image points */ 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, src_gam->cent); + rimr = icmNorm33(rimv, sci_gam->cent); /* Default setup a no compress or expand mapping of */ /* source space/image point to modified destination gamut. */ @@ -2108,8 +2008,8 @@ datao map_oh 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 = src_gam; - smp[i].dgam = dst_gam; + 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)); @@ -2143,7 +2043,7 @@ datao map_oh double tc[3] = { 0.0, 0.0, 0.0 }; for (ix = 0; ix < nmpts; ix++) { - /* Compute a rotation that brings the target point location to 50,0,0 */ + /* 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 */ @@ -2152,149 +2052,171 @@ datao map_oh } /* Figure out which neighbors of the source values to use */ - /* for the relative error & smoothing calculations. */ + /* 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++) { - int sit; - double rr; - double rrdl, rrdh; + 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",rrdl, rrdh); - +//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; - rr = sqrt(smp[ix].sv[1] * smp[ix].sv[1] + smp[ix].sv[2] * smp[ix].sv[2]); - - if (rr < 5.0) - rr = 5.0; - rr = sqrt(rr / 50.0); - - // Scale radius aprox. by cylindrical distance ?? */ - //rrdh *= rr; - - rrdl = 1.0/rrdl; - rrdh = 1.0/rrdh; - - smp[ix].nnd = 0; - - /* Until we get a minimum number of neighbors */ - for (sit = 0; smp[ix].nnd < 8 && sit < 10; sit++) { - - smp[ix].nnd = 0; - - /* Search for points within the radius */ - for (i = 0; i < nmpts; i++) { - double tt, dd, tv; - - /* Dot of neighbor color and point */ - tv = smp[i].sv[1] * smp[ix].sv[1] + smp[i].sv[2] * smp[ix].sv[2]; - - /* Ignore if of the opposote hue */ - if (tv < 0.0) - continue; - - dd = 0.0; - tt = rrdl * (smp[i].sv[0] - smp[ix].sv[0]); - dd += tt * tt; - tt = rrdh * (smp[i].sv[1] - smp[ix].sv[1]); - dd += tt * tt; - tt = rrdh * (smp[i].sv[2] - smp[ix].sv[2]); - dd += tt * tt; - - /* If we're within the filtering radius, */ - /* and not of the opposite hue */ - if (dd <= 1.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 (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); - free_nearsmth(smp, nmpts); - *npp = 0; - return NULL; - } - smp[ix].nd = nd; - smp[ix]._nnd = _nnd; + /* 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 tangent alignment 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[smp[ix].nnd].n = &smp[i]; + smp[ix].nd = nd; + smp[ix]._nnd = _nnd; + } + smp[ix].nd[smp[ix].nnd].n = &smp[i]; - /* Box filter */ -// w = 1.0; + /* Box filter */ +// w = 1.0; - /* Triangle filter */ -// w = 1.0 - dd; + /* Triangle filter */ +// w = 1.0 - dd; -// /* Cubic spline filter (default) */ - w = 1.0 - dd; - w = w * w * (3.0 - 2.0 * w); +// /* Cubic spline filter */ +// w = 1.0 - dd; +// w = w * w * (3.0 - 2.0 * w); - /* Gaussian filter */ -// w = exp(-9.0 * dd/2.0); + /* Gaussian filter (default) */ + w = exp(-9.0 * dd/2.0); - /* Sphere filter */ -// w = sqrt(1.0 - dd * dd); + /* 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; + /* Sinc^2 filter */ +// w = 3.1415926 * dd; +// if (w < 1e-9) +// w = 1e-9; +// w = sin(w)/w; +// w = w * w; - /* Save weighting */ - smp[ix].nd[smp[ix].nnd].w = w; /* Will be normalized to sum to 1.0 */ + 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); +// /* Sphere filter for depth */ +// w = sqrt(1.0 - dd * dd); - /* Cubic spline filter for depth (default) */ -// w = 1.0 - dd; -// w = w * w * (3.0 - 2.0 * w); + /* Cubic spline filter for depth */ +// w = 1.0 - dd; +// w = w * w * (3.0 - 2.0 * w); -// /* Gaussian filter for depth */ -// w = exp(-9.0 * dd/2.0); + /* Gaussian filter for depth (default) */ + w = exp(-9.0 * dd/2.0); - /* Save weighting */ - smp[ix].nd[smp[ix].nnd].rw = w; /* Won't be normalized */ + 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++; - } + smp[ix].nnd++; } - /* Increase radius in case we haven't found enough neighbors */ - rrdl /= 1.5; - rrdh /= 1.5; } - -//if (smp[ix].nnd < 8) printf("~1 point %d has %d neighbors\n",ix,smp[ix].nnd); - if (smp[ix].nnd < minnd) minnd = smp[ix].nnd; avgnd += (double)smp[ix].nnd; -//printf("~1 total of %d dir neigbours after try %d\n",smp[ix].nnd, sit); +//printf("~1 total of %d dir neigbours\n\n",smp[ix].nnd); + } avgnd /= (double)nmpts; @@ -2311,44 +2233,16 @@ datao map_oh for (j = 0; j < smp[i].nnd; j++) { smp[i].nd[j].w /= tw; } - } - } - -#ifdef SHOW_NEIGB - { - 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 magenta[3] = { 1.0, 0.0, 1.0 }; - double pp[3]; - for (i = 0; i < nmpts; i++) { - - if ((wrl = new_vrml("neigb", 1, vrml_lab)) == NULL) - error("New %s failed for '%s%s'",vrml_format(),"neigb",vrml_ext()); - for (j = 0; j < smp[i].nnd; j++) { - if (smp[i].nd[j].n == &smp[i]) - continue; - wrl->add_col_vertex(wrl, 0, smp[i].sv, yellow); - wrl->add_col_vertex(wrl, 0, smp[i].nd[j].n->sv, yellow); - } - wrl->make_lines(wrl, 0, 2); - - wrl->add_marker(wrl, smp[i].sv, red, 0.5); - - wrl->del(wrl); - printf("Waiting for input after writing 'neigb%s' for point %d:\n",vrml_ext(),i); - getchar(); } } -#endif /* SHOW_NEIGB */ #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++) { @@ -2364,7 +2258,7 @@ datao map_oh } 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); + 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); @@ -2421,8 +2315,6 @@ datao map_oh smp[i].sv[1] = smp[i].drv[1]; smp[i].sv[2] = smp[i].drv[2]; } - opts.wngam = smp[i].dgam; /* Nearest to dgam */ - opts.wn = smp[i].sv; /* minimize optfunc1 sv -> dgam */ /* Convert our start value from 3D to 2D for speed. */ icmMul3By3x4(iv, smp[i].m2d, smp[i].dv); @@ -2458,10 +2350,10 @@ datao map_oh nv[1] = iv[1] = iv[2]; powell(NULL, 2, nv, s, 0.01, 1000, optfunc1, (void *)(&opts), NULL, NULL); #endif - if (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); + 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; @@ -2499,12 +2391,14 @@ datao map_oh smp[i].aodv[2] = smp[i].drv[2]; } } + if (verb) { + printf("."); fflush(stdout); + } } - VA(("Locating weighted mapping vectors without smoothing\n")); - + VA(("Locating weighted mapping vectors without smoothing:\n")); /* Second pass to locate the optimized overall weighted point nrdv[], */ - /* which is a balance of absolute error, radial error, depth room weighting */ + /* not counting relative error. */ { double s[2] = { 20.0, 20.0 }; /* 2D search area */ double iv[3]; /* Initial start value */ @@ -2561,10 +2455,10 @@ datao map_oh nv[1] = iv[1] = iv[2]; powell(NULL, 2, nv, s, 0.01, 1000, optfunc2, (void *)(&opts), NULL, NULL); #endif - if (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); + 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; @@ -2579,90 +2473,37 @@ datao map_oh /* Remap it to the destinaton gamut surface */ smp[i].dgam->radial(smp[i].dgam, tp, tp); - icmCpy3(smp[i].dv, tp); /* Default current solution */ 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]); } - } - - /* Make sure the input and output ranges encompas the points */ - for (i = 0; i < nmpts; i++) { - for (j = 0; j < 3; j++) { - if (smp[i]._sv[j] < map_il[j]) - map_il[j] = smp[i]._sv[j];; - if (smp[i]._sv[j] > map_ih[j]) - map_ih[j] = smp[i]._sv[j]; - - if (smp[i].sv[j] < map_il[j]) - map_il[j] = smp[i].sv[j];; - if (smp[i].sv[j] > map_ih[j]) - map_ih[j] = smp[i].sv[j]; - - if (smp[i].dv[j] < map_ol[j]) - map_ol[j] = smp[i].dv[j];; - if (smp[i].dv[j] > map_oh[j]) - map_oh[j] = smp[i].dv[j]; - } - } - -#ifdef NEVER - if (verb) { - printf("Input bounding box:\n"); - printf(" %f -> %f, %f -> %f, %f -> %f\n", - map_il[0], map_ih[0], map_il[1], map_ih[1], map_il[2], map_ih[2]); - } -#endif - - /* Expand the bounding box by gexp so that our surface grid points */ - /* establish the extrapolation behaviour. Ensure that boundary */ - /* lands on the new grid though. */ - { - double scale; - - dmapres = (int)(((mapres-1) - (mapres-1)/gexp)/2.0 + 0.5); - if (dmapres < 1) - dmapres = 1; - - scale = (double)(mapres-1-dmapres)/(double)(mapres-1 - 2 * dmapres); - - for (j = 0; j < 3; j++) { - double low, high; - high = map_ih[j]; - low = map_il[j]; - map_ih[j] = (scale * (high - low)) + low; - map_il[j] = (scale * (low - high)) + high; - } -#ifdef NEVER if (verb) { - printf("After scaling up by %f, input bounding box:\n",scale); - printf(" %f -> %f, %f -> %f, %f -> %f\n", - map_il[0], map_ih[0], map_il[1], map_ih[1], map_il[2], map_ih[2]); + printf("."); fflush(stdout); } -#endif - - /* Values for grid surface points */ - hmapres = (mapres+1)/2; - hdmapres = (dmapres+1)/2; } -#if RSPLPASSES > 0 || VECADJPASSES > 0 - - VA(("Computing fine tuning correction direction:\n")); - - /* We need inward pointing correction vectors to be able */ - /* to do clipping and fine tuning. We create a shrunken */ - /* version of the dst_gamut and a mapping based on the */ - /* weighted minimum absolute error metric, and then */ - /* create a rspl to represent that mapping. */ +#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. */ - /* This sort of clipping direction helps preserve the */ - /* mapping shape (hence smoothness), while minimizing the */ - /* loss of saturation and change in dest. mapping location. */ +#if VECADJPASSES > 0 || RSPLPASSES > 0 + /* We will need inward pointing correction vectors */ { - gamut *shgam; /* Shrunken dst_gam */ + gamut *shgam; /* Shrunken di_gam */ double cusps[6][3]; double wp[3], bp[3], kp[3]; double p[3], p2[3], rad; @@ -2682,13 +2523,14 @@ datao map_oh double avgdev[MXDO]; /* Create a gamut that is a shrunk version of the destination */ - if ((shgam = new_gamut(dst_gam->getsres(dst_gam), dst_gam->getisjab(dst_gam), - dst_gam->getisrast(dst_gam))) == NULL) { + + if ((shgam = new_gamut(di_gam->getsres(di_gam), di_gam->getisjab(di_gam), + di_gam->getisrast(di_gam))) == NULL) { fprintf(stderr, "new_gamut failed\n"); - if (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); + 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; @@ -2700,14 +2542,14 @@ datao map_oh for (i = 0;;) { double len; - if ((i = dst_gam->getrawvert(dst_gam, p, i)) < 0) + if ((i = di_gam->getrawvert(di_gam, p, i)) < 0) break; doshrink(&opts, p, p, SHRINK); shgam->expand(shgam, p); } /* Translate cusps */ - if (dst_gam->getcusps(dst_gam, cusps) == 0) { + 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); @@ -2716,7 +2558,7 @@ datao map_oh shgam->setcusps(shgam, 2, NULL); } /* Translate white and black points */ - if (dst_gam->getwb(dst_gam, wp, bp, kp, NULL, NULL, NULL) == 0) { + 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); @@ -2726,10 +2568,10 @@ datao map_oh if ((gpnts = (cow *)malloc(nmpts * sizeof(cow))) == NULL) { fprintf(stderr,"gamut map: Malloc of near smooth points failed\n"); shgam->del(shgam); - if (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); + 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; @@ -2737,7 +2579,7 @@ datao map_oh /* Now locate the closest points on the shrunken gamut */ /* and set them up for creating a rspl */ - opts.wngam = shgam; + opts.shgam = shgam; for (i = 0; i < nmpts; i++) { /* Move all the points */ gtri *ctri = NULL; double tmp[3]; @@ -2749,10 +2591,9 @@ datao map_oh opts.pass = 0; /* Itteration pass */ opts.ix = i; /* Point to optimise */ opts.p = &smp[i]; - opts.wn = smp[i].dv; /* minimize optfunc1a dv -> shgam */ /* Convert our start value from 3D to 2D for speed. */ - icmMul3By3x4(iv, smp[i].m2d, smp[i].nrdv); + icmMul3By3x4(iv, smp[i].m2d, smp[i].dv); nv[0] = iv[0] = iv[1]; nv[1] = iv[1] = iv[2]; @@ -2778,16 +2619,16 @@ datao map_oh #ifdef DEBUG_POWELL_FAILS /* Optimise the point with debug on */ opts.debug = 1; - icmMul3By3x4(iv, smp[i].m2d, smp[i].nrdv); + 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, optfunc1a, (void *)(&opts), NULL, NULL); #endif shgam->del(shgam); /* Done with this */ - if (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); + 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; @@ -2803,11 +2644,12 @@ datao map_oh shgam->radial(shgam, tp, tp); /* Compute mapping vector from dst to shdst */ - icmSub3(smp[i].temp, tp, smp[i].nrdv); + 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].nrdv, &ctri); + 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); @@ -2815,12 +2657,13 @@ datao map_oh icmNormalize3(smp[i].temp, smp[i].temp, 1.0); /* Place it in rspl setup array */ - icmCpy3(gpnts[i].p, smp[i].nrdv); + 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; } @@ -2829,6 +2672,7 @@ datao map_oh 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; @@ -2846,7 +2690,7 @@ datao map_oh printf("###### gamut/nearsmth.c: writing diagnostic evects%s\n",vrml_ext()); if ((wrl = new_vrml("evects", doaxes, vrml_lab)) == NULL) error("new_vrml failed for '%s%s'","evects",vrml_ext()); - wrl->make_gamut_surface_2(wrl, dst_gam, 0.6, 0, cc); + 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); @@ -2854,16 +2698,16 @@ datao map_oh wrl->start_line_set(wrl, 0); for (i = 0; i < nmpts; i++) { - wrl->add_col_vertex(wrl, 0, smp[i].nrdv, red); + 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].nrdv, tmp); + icmSub3(tmp, smp[i].dv, tmp); #else /* Plot interpolated vectors */ - icmCpy3(cp.p, smp[i].nrdv); + icmCpy3(cp.p, smp[i].dv); evectmap->interp(evectmap, &cp); icmScale3(tmp, cp.v, 4.0); - icmSub3(tmp, smp[i].nrdv, tmp); + icmSub3(tmp, smp[i].dv, tmp); #endif wrl->add_col_vertex(wrl, 0, tmp, green); } @@ -2874,265 +2718,267 @@ datao map_oh shgam->del(shgam); /* Done with this */ free(gpnts); } +#endif /* VECADJPASSES > 0 || RSPLPASSES > 0 */ -#endif /* 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 */ - /* Lookup correction vectors */ - VA(("Smoothing guide vectors:\n")); - { - int pncliped = nmpts; - double delta; + 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 the source to destination neighborhood scale factors */ + /* 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 tmp[3]; - double sav[3], dav[3]; /* Average center locations */ - double sdev[3], ddev[3]; /* Average devation in each direction from center */ - double scev, dcev; /* Average spherical deviation */ + double cvec[3], clen; + double minext = 1e80; + double maxext = -1e80; /* Max weighted depth extension */ + double dext, gain; - for (j = 0; j < 3; j++) - sav[j] = dav[j] = sdev[j] = ddev[j] = 0.0; - scev = dcev = 0.0; + minext = -20.0; - /* Compute center average values */ - 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 */ - - icmScale3(tmp, np->sv, nw); - icmAdd3(sav, sav, tmp); - icmScale3(tmp, np->dv, nw); - icmAdd3(dav, dav, tmp); - } - - /* Compute average deviation in each direction */ + /* 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].w; /* Weight */ - double tt; - - icmSub3(tmp, sav, np->sv); - icmAbs3(tmp, tmp); - icmScale3(tmp, tmp, nw); - icmAdd3(sdev, sdev, tmp); - - tt = icmNorm33(sav, np->sv); - tt *= nw; - scev += tt; - - icmSub3(tmp, dav, np->dv); - icmAbs3(tmp, tmp); - icmScale3(tmp, tmp, nw); - icmAdd3(ddev, ddev, tmp); - - tt = icmNorm33(dav, np->dv); - tt *= nw; - dcev += tt; - } - -//printf("~1 %d: sdev %f %f %f, scev %f\n",i,sdev[0],sdev[1],sdev[2],scev); -//printf("~1 %d: ddev %f %f %f, dcev %f\n",i,ddev[0],ddev[1],ddev[2],dcev); + double nw = smp[i].nd[j].rw; /* Weight */ + double tmpl; - /* Try and protect against silliness */ - if (scev < 1e-3 || dcev < 1e-3) - scev = dcev = 1e-3; + icmSub3(cvec, np->tdst, np->anv); /* Vector needed to target for neighbour */ + clen = icmDot3(smp[i].evect, cvec); /* Error in this direction */ - for (j = 0; j < 3; j++) { - if (sdev[j] < 1e-3 || ddev[j] < 1e-3) { - sdev[j] = scev; - ddev[j] = dcev; - } + tmpl = nw * (clen - minext); /* Track maximum weighted extra depth */ + if (tmpl < 0.0) + tmpl = 0.0; + if (tmpl > maxext) + maxext = tmpl; } + maxext += minext; - /* Compute scale factors */ - icmDiv3(smp[i].nscale, ddev, sdev); /* Scale = ddev/sdev */ - -#ifdef NEVER -if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 - || smp[i].nscale[1] > 1.5 || smp[i].nscale[1] < 0.01 - || smp[i].nscale[2] > 1.5 || smp[i].nscale[2] < 0.01) { - printf("~1 %d: scale factors %f %f %f\n",i,smp[i].nscale[0], smp[i].nscale[1], smp[i].nscale[2]); - printf("~1 %d: from sdev %f %f %f\n",i,sdev[0], sdev[1], sdev[2]); - printf("~1 %d: from ddev %f %f %f\n",i,ddev[0], ddev[1], ddev[2]); -} -#endif /* NEVER */ + if (it == 0) + gain = 1.2; + else + gain = 0.8; - } + /* Accumulate correction with damping */ + smp[i].rext += gain * maxext; - /* Itterate smoothing until we're happy */ - for (it = 0; it < VECADJPASSES; it++) { - int ncliped = 0; - double maxclipby = 0.0; - double avgclipby = 0.0; - - /* Compute the neighbourhood smoothed anv[] from dv[] */ - for (i = 0; i < nmpts; i++) { - double sav[3], dav[3]; /* Average locations */ - double tmp[3], c1[3], c2[3]; - double rdsm; - - /* Compute average values */ - sav[0] = sav[1] = sav[2] = 0.0; - dav[0] = dav[1] = dav[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 */ - - icmScale3(tmp, np->sv, nw); /* weight for filter */ - icmAdd3(sav, sav, tmp); /* sum filtered value */ - - /* weight for filter */ - tmp[0] = nw * np->dv[0]; /* Don't itterate J */ - tmp[1] = nw * np->anv[1]; - tmp[2] = nw * np->anv[2]; - icmAdd3(dav, dav, tmp); /* sum filtered value */ - } - - /* Compute filtered value with source to dest scaling */ - icmSub3(tmp, smp[i].sv, sav); /* Vector from average to src */ - icmMul3(tmp, tmp, smp[i].nscale); /* Scale */ - icmAdd3(tmp, tmp, dav); /* average dst + vector */ + /* Error for just this point */ + icmSub3(cvec, smp[i].tdst, smp[i].anv); + clen = icmDot3(smp[i].evect, cvec); - rdsm = 1.0 - sqrt(smp[i].wt.r.dsm); /* To degree of blending with unchanged */ + /* Blend to individual correction on neutral axis */ + dext = smp[i].naxbf * smp[i].rext + (1.0 - smp[i].naxbf) * clen; - icmBlend3(tmp, tmp, smp[i].dv, rdsm); /* Less than full imprint */ - -#if VECADJPASSES > 1 - /* Clip to gamut */ - if (dc_gam->nradial(dc_gam, c1, tmp) > (1.0 + 1e-6)) { - co cp; - double cvec[3]; - - /* Lookup "shrunk gamut" cliping direction */ - icmCpy3(cp.p, tmp); - evectmap->interp(evectmap, &cp); - icmNormalize3(cvec, cp.v, 1.0); - - if (!vintersect2(dc_gam, NULL, c2, cvec, tmp)) { /* Got an intersection */ - double id; - -//printf("~1 clipped %f %f %f -> %f %f %f\n", tmp[0], tmp[1], tmp[2], c2[0], c2[1], c2[2]); - id = icmNorm33(c2, tmp); /* Dist to intersection */ - icmCpy3(tmp, c2); - - ncliped++; - if(id > maxclipby) - maxclipby = id; - avgclipby += id; - } else { -//printf("~1 rclipped %f %f %f -> %f %f %f\n", tmp[0], tmp[1], tmp[2], c1[0], c1[1], c1[2]); - icmCpy3(tmp, c1); /* Use radial clip */ - } - } -#endif - - /* Blend to un-smoothed value on neutral axis */ - icmBlend3(tmp, smp[i].dv, tmp, smp[i].naxbf); - - /* Updated value for next itteration */ - icmCpy3(smp[i].anv, tmp); - } + /* Apply integrated correction */ + icmScale3(cvec, smp[i].evect, dext); + icmAdd3(smp[i].anv, smp[i].dv, cvec); - if (ncliped > 0) - avgclipby /= (double)ncliped; + if (clen > 0.0) { /* Compression */ + if (clen > maxog) + maxog = clen; + avgog += clen; + nog++; - delta = (pncliped - ncliped)/(double)nmpts; - - if (verb) { - printf("It %d: No clip %d/%d delta %f max by %f, avg by %f\n",it,ncliped, nmpts+1, delta, maxclipby, avgclipby); + } else { /* Expansion */ + if (-clen > maxig) + maxig = -clen; + avgig += -clen; + nig++; } - pncliped = ncliped; - } - - /* 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) + 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 /* VECADJPASSES > 0 */ - -#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 */ - smp[i].dr = icmNorm33(smp[i].dv, smp[i].dgam->cent); /* Vector smoothed */ + /* 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); } -#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 RSPLPASSES > 0 + if (verb) { + double avgog = 0.0, maxog = 0.0, nog = 0.0; + double avgig = 0.0, maxig = 0.0, nig = 0.0; - VA(("Fine tuning vectors to allow for rspl smoothing:\n")); + /* 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"); if (evectmap != NULL) evectmap->del(evectmap); - if (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); + 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; } - /* Lookup correction vectors */ - VA(("Computing fine tuning target for vectors:\n")); + /* Loopkup correction vectors */ + VA(("Computing fine tuning direction for vectors:\n")); for (i = 0; i < nmpts; i++) { + co cp; double nd, id, tmp[3]; - /* If the sv and dv are within dc_gam, then this point doesn't need */ - /* to be fine tuned to make it land on the gamut surface - this point */ - /* either doesn't need gamut mapping, or is being expanded, in which */ - /* case we prioritize smoothness over exactly hitting the expansion */ - /* target */ - if (dc_gam->nradial(dc_gam, NULL, smp[i].sv) <= (1.0 + 1e-6) - && dc_gam->nradial(dc_gam, NULL, smp[i].dv) <= (1.0 + 1e-6)) { - icmCpy3(smp[i].tdst, smp[i].dv); /* Target is where we are */ - smp[i].nott = 1; - - } else { - co cp; - double evect[3]; - - /* Lookup fine tuning vector direction for current location */ - icmCpy3(cp.p, smp[i].dv); - evectmap->interp(evectmap, &cp); - icmNormalize3(evect, cp.v, 1.0); - - /* 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, 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].nott = 0; + 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; @@ -3150,8 +2996,6 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 double avgrext = 0.0; double ovlen; - VA(("it %d: Creating rspl\n",it)); - /* 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 */ @@ -3160,14 +3004,13 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 } 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, - map_il, map_ih, gres, map_ol, map_oh, mapsmooth, avgdev, NULL); - - VA(("it %d: Evaluate mapping\n",it)); + 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. */ @@ -3179,20 +3022,12 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 icmCpy3(cp.p, smp[i]._sv); map->interp(map, &cp); icmCpy3(smp[i].temp, cp.v); - - /* Lookup fine tuning vector direction for that value. */ - /* (evect[] is then used in the local correction loop below) */ - icmCpy3(cp.p, smp[i].temp); - evectmap->interp(evectmap, &cp); - icmNormalize3(smp[i].evect, cp.v, 1.0); - + /* 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); } - VA(("it %d: Compute correction vectors\n",it)); - /* Compute local correction */ for (i = 0; i < nmpts; i++) { double minext = 1e80; @@ -3202,8 +3037,6 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 double tt; double cgain, xgain; /* This itters compression, expansion gain */ double gain; /* Gain used */ - co cp; - double evect[3]; /* See what the worst case is in the local area, and */ /* aim to lower the whole local area by enough to */ @@ -3258,7 +3091,6 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 /* Keep stats of this point */ clen = smp[i].clen; - if (clen > 0.0) { if (clen > maxog) maxog = clen; @@ -3283,32 +3115,20 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 gpnts[i].w = 1.0; } - if ((it+1) < RSPLPASSES || !surfpnts) - map->del(map); /* Not the last pass, or not doing grid surface points */ - else - lastmap = map; /* Let grid surface creation use this. */ - - VA(("it %d: Compute correction rspl\n",it)); - - /* Create rspl of corrections */ 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, - map_il, map_ih, gres, map_ol, map_oh, 1.0, avgdev, NULL); - - VA(("it %d: Apply corrections\n",it)); + 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; - if (smp[i].nott) /* Don't alter points within the gamut */ - continue; - icmCpy3(cp.p, smp[i].dv); map->interp(map, &cp); #ifdef RSPLUSEPOW @@ -3323,7 +3143,7 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 /* Apply accumulated offset */ icmAdd3(smp[i].anv, smp[i].dv, cp.v); } - map->del(map); /* Not the last pass, or not doing grid surface points */ + 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", @@ -3354,15 +3174,18 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 } #endif /* RSPLPASSES > 0 */ -#endif /* !DIAG_POINTS */ +#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, mapres); + create_influence_plot(smp, nmpts); #endif - VA(("Restoring non cusp-rotated source points:\n")); + VB(("Final guide points:\n")); /* Restore the actual non cusp rotated source point */ for (i = 0; i < nmpts; i++) { @@ -3413,7 +3236,7 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 /* Compute actual depth of ray into destination gamut */ /* to determine if this is expansion or contraction. */ - if (dst_gam->vector_isect(dst_gam, smp[i].sv, smp[i].dv, + 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 */ @@ -3426,7 +3249,7 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 /* 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 (dc_gam->getwb(dc_gam, NULL, NULL, NULL, wp, dst_kbp ? NULL : bp, dst_kbp ? bp : NULL) == 0) { + 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) { double nalev[3]; @@ -3472,7 +3295,7 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 } #ifdef VERB else { - printf("dc_gam->getwb failed\n"); + printf("d_gam->getwb failed\n"); } #endif @@ -3621,345 +3444,54 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 #ifdef SUBVEC_SMOOTHING VB(("Smoothing sub-surface guide points:\n")); - { - double maxmv = 0.0, avgmv = 0.0, acount = 0.0; - /* Smooth the sub-surface mapping points */ - for (i = 0; i < nmpts; i++) { - double sav[3], dav[3]; /* Average locations */ - double scr, dcr; /* Cylindrical radius */ - double scf; /* Scale factor */ - double tmp[3], de; - - if (smp[i].vflag == 0) - continue; /* Sub value not valid */ - - /* Compute average values */ - sav[0] = sav[1] = sav[2] = 0.0; - dav[0] = dav[1] = dav[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 */ + /* 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]; - icmScale3(tmp, np->sv2, nw); /* weight for filter */ - icmAdd3(sav, sav, tmp); /* sum filtered value */ + 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, np->dv2, nw); /* weight for filter */ - icmAdd3(dav, dav, tmp); /* sum filtered value */ + icmScale3(tmp, tmp, nw); /* weight for filter */ + icmAdd3(fdv2, fdv2, tmp); /* sum filtered value */ + tw += nw; } - - /* We want to transfer the relative location (i.e. detail) from */ - /* the source to destination, but we need to scale the features */ - /* appropriately for the mapping. */ - scr = sqrt(sav[1] * sav[1] + sav[2] * sav[2]); - dcr = sqrt(dav[1] * dav[1] + dav[2] * dav[2]); - scf = dcr/scr; - - /* Compute filtered value */ - icmSub3(tmp, smp[i].sv2, sav); /* Vector from average to src */ - tmp[1] *= scf; /* Scale */ - tmp[2] *= scf; /* Scale */ - icmAdd3(tmp, tmp, dav); /* average dst + vector */ - - de = icmNorm33(smp[i].dv2, tmp); - icmCpy3(smp[i].dv2, tmp); - - if (de > maxmv) - maxmv = de; - avgmv += de; - acount++; - - VB(("Smthd Src %d = %f %f %f\n",i,smp[i].sv2[0],smp2[i].sv[1],smp2[i].sv2[2])); - VB(("Smthd Dst %d = %f %f %f\n",i,smp[i].dv2[0],smp2[i].dv[1],smp2[i].dv2[2])); } - if (acount > 0) - avgmv /= acount; - - if (verb) - printf("Sub-surface smoothing changed by max %f, average %f\n",maxmv, avgmv); + 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")); -#ifdef PLOT_SMOOTHING_CHANGE - /* Plot change in destination point of un-smoothed to smoothed */ - { - vrml *wrl = NULL; - int doaxes = 0; - -#ifdef PLOT_AXES - doaxes = 1; -#endif - wrl = new_vrml("dst_smvec", doaxes, vrml_lab); - - /* Start of guide vector plot */ - wrl->start_line_set(wrl, 0); - - for (i = 0; i < nmpts; i++) { - double red[3] = { 1.0, 0.0, 0.0 }; - double green[3] = { 0.0, 1.0, 0.0 }; - - wrl->add_col_vertex(wrl, 0, smp[i].nrdv, red); - wrl->add_col_vertex(wrl, 0, smp[i].dv, green); - } - wrl->make_lines(wrl, 0, 2); /* Change vectors */ - -#ifndef NEVER - /* Plot un-smoothed src to dst mappings */ - wrl->start_line_set(wrl, 0); - - for (i = 0; i < nmpts; i++) { - double lblue[3] = { 0.4, 0.4, 0.8 }; - double magenta[3] = { 0.8, 0.4, 0.8 }; - - wrl->add_col_vertex(wrl, 0, smp[i].sv, lblue); - wrl->add_col_vertex(wrl, 0, smp[i].nrdv, magenta); - } - wrl->make_lines(wrl, 0, 2); /* Change vectors */ -#endif - -#ifdef NEVER - /* Plot index numbers */ - for (i = 0; i < nmpts; i++) { - double cream[3] = { 0.7, 0.7, 0.5 }; - char buf[100]; - sprintf(buf, "%d", i); - wrl->add_text(wrl, buf, smp[i].dv, cream, 0.5); - } -#endif /* NEVER */ - - /* Write transparent destination space gamut surface */ - dc_gam->write_to_vrml(dc_gam, wrl, 0.5, 0); - - /* Write file */ - wrl->del(wrl); - } -#endif /* PLOT_SMOOTHING_CHANGE */ - - /* If grid surface points are requested */ - if (surfpnts) { - DCOUNT(gc, 3, 3, 0, 0, hmapres); - double cent[3]; - - VB(("Adding grid surface points:\n")); - - /* If rspl smoothing didn't leave us a map */ - if (lastmap == NULL) { - - cow *gpnts = NULL; /* Mapping points to create 3D -> 3D mapping */ - int gres[MXDI]; - double avgdev[MXDO]; - - VB(("Creating rspl map for grid surface points\n",it)); - - if ((gpnts = (cow *)malloc(nmpts * sizeof(cow))) == NULL) { - fprintf(stderr,"gamut map: Malloc of near smooth points failed\n"); - if (evectmap != NULL) - evectmap->del(evectmap); - if (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); - free_nearsmth(smp, nmpts); - *npp = 0; - return NULL; - } - - /* Setup the rspl guide points for creating rspl */ - for (i = 0; i < nmpts; i++) { - icmCpy3(gpnts[i].p, smp[i].sv); - icmCpy3(gpnts[i].v, smp[i].dv); - gpnts[i].w = 1.0; - } - - for (j = 0; j < 3; j++) { /* Set resolution for all axes */ - gres[j] = mapres; /* Full resolution */ - avgdev[j] = GAMMAP_RSPLAVGDEV; - } - lastmap = new_rspl(RSPL_NOFLAGS, 3, 3); /* Allocate 3D -> 3D */ - lastmap->fit_rspl_w(lastmap, GAMMAP_RSPLFLAGS, gpnts, nmpts, - map_il, map_ih, gres, map_ol, map_oh, mapsmooth, avgdev, NULL); - free(gpnts); - } - - sc_gam->getcent(dc_gam, cent); - - DC_INIT(gc); - for (;;) { - /* If point is in the outer two layers of grid */ - if ( gc[0] == 0 || gc[0] == hdmapres - || gc[0] == (hmapres-1) || gc[0] == (hmapres-1-hdmapres) - || gc[1] == 0 || gc[1] == hdmapres - || gc[1] == (hmapres-1) || gc[1] == (hmapres-1-hdmapres) - || gc[2] == 0 || gc[2] == hdmapres - || gc[2] == (hmapres-1) || gc[2] == (hmapres-1-hdmapres)) - - /* Only points around gamut, not on top or underneath */ -/* - if ( gc[1] == 0 || gc[1] == hdmapres - || gc[1] == (hmapres-1) || gc[1] == (hmapres-1-hdmapres) - || gc[2] == 0 || gc[2] == hdmapres - || gc[2] == (hmapres-1) || gc[2] == (hmapres-1-hdmapres)) -*/ - { - double grid2gamut, gamut2cent, ww; - co cp; - - if (nmpts >= mxnmpts) { - warning("nearsmth ran out of space for points"); - break; - } - smp[nmpts].uflag = 1; - - /* Source location */ - for (j = 0; j < 3; j++) - smp[nmpts].sv[j] = map_il[j] + gc[j]/(hmapres-1.0) * (map_ih[j] - map_il[j]); - - /* If this point is within source gamut, skip it */ - if (sc_gam->nradial(sc_gam, NULL, smp[nmpts].sv) <= (1.0 + 1e-6)) { -//printf("~1 point %d %d %d = %f %f %f is inside source gamut\n", gc[0], gc[1], gc[2], smp[nmpts].sv[0], smp[nmpts].sv[1], smp[nmpts].sv[2]); - goto next_point; - } -#ifdef NEVER - /* Clip the point to the closest location on the source */ - /* colorspace gamut. */ - sc_gam->nearest(sc_gam, cp.p, smp[nmpts].sv); -#else - /* Map grid point to weighted nearest on source space gamut */ - { - double ta[3] = { 50.0, 0.0, 0.0 }; - double tc[3] = { 0.0, 0.0, 0.0 }; - double s[2] = { 20.0, 20.0 }; /* 2D search area */ - double nv[2]; /* 2D New value */ - double tp[3]; /* Resultint value */ - double ne; /* New error */ - int notrials = NO_TRIALS; - double bnv[3]; /* Best 3d value */ - double brv; /* Best return value */ - int trial; - double mv; - - /* Determine the parameter weighting at this location */ - opts.pass = 0; /* Itteration pass */ - opts.ix = nmpts; - opts.p = &smp[nmpts]; - opts.wngam = sc_gam; /* Optimise to source colorspace gamut */ - opts.wn = smp[nmpts].sv; /* minimize optfunc1a sv -> sc_gam */ - - /* Compute weights at this point */ - interp_xweights(sc_gam, &smp[nmpts].wt, smp[nmpts].sv, opts.xwh, &opts, 0); - - /* Initial starting point */ - sc_gam->nearest(sc_gam, bnv, smp[nmpts].sv); - - /* 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 */ - - /* Setup the 3D -> 2D tangent conversion and inverse for our start point */ - icmVecRotMat(smp[nmpts].m2d, bnv, sc_gam->cent, ta, tc); - icmVecRotMat(smp[nmpts].m3d, ta, tc, bnv, sc_gam->cent); - - /* Convert our start value from 3D to 2D for speed. */ - icmMul3By3x4(tp, smp[nmpts].m2d, bnv); - nv[0] = tp[1]; - nv[1] = tp[2]; - - if (trial >= 2) { - /* Use random offset to avoid local minima */ - nv[0] += d_rand(-20.0, 20.0); - nv[1] += d_rand(-20.0, 20.0); - } - - /* Optimise the point */ - if (powell(&rv, 2, nv, s, 0.01, 1000, optfunc1a, (void *)(&opts), NULL, NULL) == 0 - && rv < brv) { - brv = rv; -//printf("~1 point %d, trial %d, new best %f\n",i,trial,rv); - - /* Convert best result 2D -> 3D */ - tp[2] = nv[1]; - tp[1] = nv[0]; - tp[0] = 50.0; - icmMul3By3x4(tp, smp[nmpts].m3d, tp); - - /* Remap it to the source gamut surface */ - sc_gam->radial(sc_gam, bnv, tp); - } -//else printf("~1 powell failed with rv = %f\n",rv); - } - if (brv == 1e38) { /* We failed to get a result */ - fprintf(stderr, "multiple powells failed to get a result (4)\n"); - sc_gam->nearest(sc_gam, cp.p, smp[nmpts].sv); - - } else { - icmCpy3(cp.p, bnv); - } - } -#endif /* NEVER */ - -//printf("~1 grid %f %f %f -> src %f %f %f\n", smp[nmpts].sv[0], smp[nmpts].sv[1], smp[nmpts].sv[2], cp.p[0], cp.p[1], cp.p[2]); - - /* Then lookup the gamut mapped value */ - lastmap->interp(lastmap, &cp); - -//printf("~1 src %f %f %f -> dst %f %f %f\n", cp.p[0], cp.p[1], cp.p[2], cp.v[0], cp.v[1], cp.v[2]); - - for (j = 0; j < 3; j++) - smp[nmpts].dv[j] = cp.v[j]; - - /* Compute the distance of the grid surface point to the to the */ - /* source colorspace gamut, as well as the distance from there */ - /* to the gamut center point. */ - for (grid2gamut = gamut2cent = 0.0, j = 0; j < 3; j++) { - double tt; - tt = smp[nmpts].dv[j] - cp.p[j]; - grid2gamut += tt * tt; - tt = cp.p[j] - cent[j]; - gamut2cent += tt * tt; - } - grid2gamut = sqrt(grid2gamut); - gamut2cent = sqrt(gamut2cent); - if (gamut2cent < 0.1) - gamut2cent = 0.1; - - /* Make the weighting inversely related to distance, */ - /* to reduce influence on in gamut mapping shape, */ - /* while retaining some influence at the edge of the */ - /* grid. */ - ww = grid2gamut / gamut2cent; - if (ww > 1.0) - ww = 1.0; - - /* A low weight seems to be enough ? */ - /* The lower the better in terms of geting best hull mapping fidelity */ - smp[nmpts++].w1 = 0.1 * ww; - } - next_point:; - DC_INC(gc); - if (DC_DONE(gc)) - break; - } - *npp = nmpts; /* Update returned number of points */ - - lastmap->del(lastmap); - } - if (evectmap != NULL) evectmap->del(evectmap); #ifndef PLOT_DIGAM - if (src_gam != sc_gam) - src_gam->del(src_gam); - if (dst_gam != src_gam && dst_gam != dc_gam) - dst_gam->del(dst_gam); + 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; - smp[i].dcgam = NULL; } #else /* !PLOT_DIGAM */ warning("!!!!! PLOT_DIGAM defined !!!!!"); @@ -3972,7 +3504,6 @@ if (smp[i].nscale[0] > 1.5 || smp[i].nscale[0] < 0.01 void free_nearsmth(nearsmth *smp, int nmpts) { int i; - /* Free contents that have been used */ for (i = 0; i < nmpts; i++) { if (smp[i].nd != NULL) free(smp[i].nd); @@ -3987,7 +3518,7 @@ void free_nearsmth(nearsmth *smp, int nmpts) { /* 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, int mapres) { +static void create_influence_plot(nearsmth *smp, int nmpts) { int i, j, k; gamut *gam; int src = 0; /* 1 = src, 0 = dst gamuts */ @@ -4061,7 +3592,7 @@ static void create_influence_plot(nearsmth *smp, int nmpts, int mapres) { /* Create the diagnostic color rspl */ for (j = 0; j < 3; j++) { /* Set resolution for all axes */ - gres[j] = mapres; + gres[j] = smp->mapres; avgdev[j] = 0.001; } swdiag = new_rspl(RSPL_NOFLAGS, 3, 3); /* Allocate 3D -> 3D */ diff --git a/gamut/nearsmth.h b/gamut/nearsmth.h index 82afab3..f826357 100644 --- a/gamut/nearsmth.h +++ b/gamut/nearsmth.h @@ -126,14 +126,10 @@ typedef struct { struct { double rdl; /* Direction smoothing radius L* dir. (delta E radius at src point)*/ double rdh; /* Direction smoothing radius H* (delta E radius at src point)*/ - - double dsm; /* Degree of smoothing (non-linear response) */ } r; - /* Depth room weighting. */ - /* Weighing to give to minimizing depth ratio by mapping to/from adequate dest/src depth. */ - /* The idea is to compromize luminance and/or hue to allow room for */ - /* preserving saturation distinction under heavy compression. */ + /* depth weighting */ + /* Weighing to give to minimizing depth ratio by mapping to/from adequate dest/src depth */ struct { double co; /* Overall compression weighting */ double xo; /* Overall expansion weighting */ @@ -170,7 +166,6 @@ typedef struct { struct _nearsmth { /* Public: */ - int uflag; /* Use flag, 0 = normal, 1 = grid surface point */ int gflag; /* Gamut direction flag. 0 = not determinable, 1 = comp., 2 = exp. */ int vflag; /* Vector direction flag. 0 = not determinable, 1 = comp., 2 = exp. */ /* sv2, dv2, sd3 etc. are valid if vflag != 0 */ @@ -183,13 +178,12 @@ struct _nearsmth { double dv[3]; /* Output destination value */ double dr; /* Output destination value radius from center */ double div[3]; /* gam[cx]pf moderated dv[] value */ - double w1; /* guide point weight */ /* Gamut sub-surface mapping guide point (knee shape controlled by gamcknf & gamxknf) */ double sv2[3]; /* Sub-surface source value */ double dv2[3]; /* Sub-surface knee'd adjusted destination value */ double div2[3]; /* gam[cx]pf moderated dv2[] value */ - double w2; /* Sub-surface weight (set in nearsmth) */ + double w2; /* Sub-surface weight (fixed in nearsmth) */ double sd3[3]; /* Deep sub-surface source & destination value */ double w3; /* Deep sub-surface weight */ @@ -207,43 +201,39 @@ struct _nearsmth { double _sr; /* Original source radius */ double naxbf; /* Blend factor that goes to 0.0 at white & black points. */ double aodv[3]; /* Absolute error optimized destination value */ - double nrdv[3]; /* No relative smoothed destination value */ + double nrdv[3]; /* No relative weight optimized destination value */ double anv[3]; /* Average neighborhood target point (relative target) */ double tdst[3]; /* Target destination on gamut */ - int nott; /* NZ if not a point that needs to land on gamut */ double evect[3]; /* Accumulated extension vector direction */ double clen; /* Current correction length needed */ double coff[3]; /* Correction offset */ double rext; double temp[3]; /* General temporary */ - gamut *sgam; /* Source gamut src_gam = intersection of src and img gamut gamut */ - gamut *dgam; /* Intersected destination gamut dst_gam */ - gamut *dcgam; /* Destination Colorspace gamut dc_gam */ + gamut *sgam; /* Source gamut sci_gam = intersection of src and img gamut gamut */ + gamut *dgam; /* Destination gamut di_gam */ int nnd, _nnd; /* No & size of direction neighbour list */ neighb *nd; /* Allocated list of neighbours */ - double nscale[3]; /* Neighborhood scale change from sv to dv */ double dcratio; /* Depth compression ratio */ double dxratio; /* Depth expansion ratio */ + int mapres; /* Target grid res for 3D RSPL */ + int debug; - double dbgv[3]; /* Error components va, vr, vd on last itteration */ + double dbgv[4]; /* Error components va, vr, vl, vd on last itteration */ }; typedef struct _nearsmth nearsmth; /* Return the upper bound on the number of points that will be generated */ int near_smooth_np( - gamut **pp_gam, /* Return gamut that was used for points */ gamut *sc_gam, /* Source colorspace gamut */ gamut *s_gam, /* Source image gamut (== sc_gam if none) */ gamut *d_gam, /* Destination colorspace gamut */ - double xvra, /* Extra vertex ratio */ - int gmult, /* Guide point multiplier, typically 4 */ - int surfgres /* surface grid point resolution, 0 for none */ + double xvra /* Extra vertex ratio */ ); /* Return a list of points. Call free_nearsmth() after use */ @@ -263,13 +253,11 @@ nearsmth *near_smooth( 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 vertex ratio */ - int mapres, /* Target grid res for 3D RSPL, (allowing for gexp) */ + int mapres, /* Target grid res for 3D RSPL */ double mapsmooth, /* Target smoothing for 3D RSPL */ - double gexp, /* Total grid expansion ratio, none = 1.0 */ - int surfgres, /* Surface grid point resolution, 0 for none */ - datai map_il, /* Return input range */ + datai map_il, /* Preliminary rspl input range */ datai map_ih, - datao map_ol, /* Return output range */ + datao map_ol, /* Preliminary rspl output range */ datao map_oh ); @@ -279,21 +267,13 @@ void free_nearsmth(nearsmth *smp, int npp); /* Expand the compact form of weights into the explicit form. */ int expand_weights(gammapweights out[14], gammapweights *in); -/* Blend two expanded groups of individual weights into one */ +/* Blend a two expanded groups of individual weights into one */ void near_xwblend( gammapweights *dst, gammapweights *src1, double wgt1, gammapweights *src2, double wgt2 ); -/* Blend three expanded groups of individual weights into one */ -void near_xwblend3( -gammapweights *dst, -gammapweights *src1, double wgt1, -gammapweights *src2, double wgt2, -gammapweights *src3, double wgt3 -); - /* Tweak weights acording to extra cmy cusp flags or rel override */ void tweak_weights(gammapweights out[14], int dst_cmymap, int rel_oride); diff --git a/gamut/smthtest.c b/gamut/smthtest.c index b794ed0..c9214d3 100644 --- a/gamut/smthtest.c +++ b/gamut/smthtest.c @@ -47,44 +47,37 @@ gammapweights weights[] = { { /* Cusp alignment control */ { 0.1, /* Cusp luminance alignment weighting 0 = none, 1 = full */ - 0.0, /* Cusp chroma alignment weighting 0 = none, 1 = full */ - 0.3 /* Cusp hue alignment weighting 0 = none, 1 = full */ + 0.1, /* Cusp chroma alignment weighting 0 = none, 1 = full */ + 0.2 /* 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 (currently broken - need to fix) */ + { /* Radial weighting */ 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.8, /* Hue dominance vs l+c, 0 - 1 */ + 0.5, /* 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 */ - 0.93, /* Black l dominance vs, c, 0 - 1 */ + 0.9, /* Light l dominance vs, c, 0 - 1 */ + 0.9, /* Medium l dominance vs, c, 0 - 1 */ + 0.9, /* Dark l dominance vs, c, 0 - 1 */ - 0.4, /* White l blend start radius, 0 - 1, at white = 0 */ - 0.7, /* Black l blend power, linear = 1.0, enhance < 1.0 */ - - 1.5, /* L error extra power with size, none = 1.0 */ - 10.0 /* L error extra xover threshold in DE */ + 0.5, /* l/c dominance breakpoint, 0 - 1 */ + 0.0, /* l dominance exageration, 0+ */ + 0.0 /* c dominance exageration, 0+ */ }, { /* Relative vector smoothing */ - 20.0, 30.0, /* Relative Smoothing radius L* H* */ - 0.9 /* Degree of smoothing */ + 30.0, 20.0 /* Relative Smoothing radius L* H* */ }, { /* Weighting of excessive compression error, which is */ /* the src->dst vector length over the available dst depth. */ /* The depth is half the distance to the intersection of the */ /* vector to the other side of the gamut. (doesn't get triggered much ?) */ - 5.0, /* Compression depth weight */ - 5.0 /* Expansion depth weight */ - }, - { - 0.0 /* Fine tuning expansion weight, 0 - 1 */ + 100.0, /* Compression depth weight */ + 100.0 /* Expansion depth weight */ } } }; @@ -217,7 +210,7 @@ main(int argc, char *argv[]) { /* Create the near point mapping */ nsm = near_smooth(verb, &nnsm, gin, gin, gout, 0, 0, NULL, xweights, - 0.1, 0.1, 1, 1, 2.0, 19, 2.0, 1.20, 5, il, ih, ol, oh); + 0.1, 0.1, 1, 1, 2.0, 17, 10.0, il, ih, ol, oh); if (nsm == NULL) error("Creating smoothed near points failed"); |