summaryrefslogtreecommitdiff
path: root/gamut
diff options
context:
space:
mode:
Diffstat (limited to 'gamut')
-rw-r--r--gamut/gammap.c667
-rw-r--r--gamut/gamut.c682
-rw-r--r--gamut/gamut.h25
-rw-r--r--gamut/maptest.c3
-rw-r--r--gamut/nearsmth.c1803
-rw-r--r--gamut/nearsmth.h48
-rw-r--r--gamut/smthtest.c35
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");
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-14 19:22:39 +0000