summaryrefslogtreecommitdiff
path: root/gamut/gammap.c
diff options
context:
space:
mode:
authorJörg Frings-Fürst <debian@jff-webhosting.net>2014-09-01 13:56:46 +0200
committerJörg Frings-Fürst <debian@jff-webhosting.net>2014-09-01 13:56:46 +0200
commit22f703cab05b7cd368f4de9e03991b7664dc5022 (patch)
tree6f4d50beaa42328e24b1c6b56b6ec059e4ef21a5 /gamut/gammap.c
Initial import of argyll version 1.5.1-8debian/1.5.1-8
Diffstat (limited to 'gamut/gammap.c')
-rw-r--r--gamut/gammap.c2552
1 files changed, 2552 insertions, 0 deletions
diff --git a/gamut/gammap.c b/gamut/gammap.c
new file mode 100644
index 0000000..a9bef1e
--- /dev/null
+++ b/gamut/gammap.c
@@ -0,0 +1,2552 @@
+
+/*
+ * Argyll Gamut Mapping Library
+ *
+ * Author: Graeme W. Gill
+ * Date: 1/10/00
+ * Version: 2.00
+ *
+ * Copyright 2000 - 2006 Graeme W. Gill
+ * All rights reserved.
+ *
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ *
+ * For a discussion of gamut mapping strategy used,
+ * see gammap.txt
+ */
+
+/*
+ * TTBD:
+ * Improve error handling.
+ *
+ * There is a general expectation (especially in comparing products)
+ * that the profile colorimetric intent be not strictly minimum delta E,
+ * but that it correct neutral axis, luminence range and keep hue
+ * proportionality. Ideally there should be an intent that matches
+ * this, that can be selected for the colorimetric table (or perhaps be default).
+ *
+ * It might be good to offer the black mapping method as an option (gmm_BPmap),
+ * as well as offering different profile (xicc/xlut.c) black point options
+ * (neutral, K hue max density, CMY max density, any max density).
+ *
+ * The gamut mapping code here and the near smooth code don't actually mesh
+ * very well. For instance, the black point bend approach in < V1.3.4
+ * means that the dest gamut isn't actually contained within the source,
+ * messing up the guide vector mappings. Even if this is fixed, the
+ * actual neutral aim point within nearsmooth is Jab 0,0, while
+ * the mapping in gammap is from the source neutral to the chosen
+ * ??????
+ */
+
+
+
+#define VERBOSE /* [Def] Print out extra interesting information when verbose is set */
+#undef PLOT_DIAG_WRL /* [Und] Always plot "gammap.wrl" */
+
+ /* What do display when user requests disgnostic .wrl */
+#define PLOT_SRC_GMT /* [Def] Plot the source surface to "gammap.wrl" as well */
+#define PLOT_DST_GMT /* [Def] Plot the dest surface to "gammap.wrl" as well */
+#undef PLOT_SRC_CUSPS /* [Und] Plot the source surface cusps to "gammap.wrl" as well */
+#undef PLOT_DST_CUSPS /* [Und] Plot the dest surface cusps to "gammap.wrl" as well */
+#undef PLOT_TRANSSRC_CUSPS /* [Und] Plot the gamut mapped source surface cusps to "gammap.wrl" */
+#undef PLOT_AXES /* [Und] Plot the axes to "gammap.wrl" as well */
+#undef SHOW_VECTOR_INDEXES /* [Und] Show the mapping vector index numbers */
+#define SHOW_MAP_VECTORS /* [Def] Show the mapping vectors */
+#undef SHOW_SUB_SURF /* [Und] Show the sub-surface mapping vector */
+#undef SHOW_CUSPMAP /* [Und] Show the cusp mapped vectors rather than final vectors */
+#undef SHOW_ACTUAL_VECTORS /* [Und?] Show how the source vectors actually map thought xform */
+#undef SHOW_ACTUAL_VEC_DIFF /* [Und] Show how the difference between guide and actual vectors */
+
+#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 the 3D compression knees */
+#undef CHECK_NEARMAP /* [Und] Check how accurately near map vectors are represented by rspl */
+
+#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 */
+
+#undef PLOT_DIGAM /* [Und] Rather than DST_GMT - don't free it (#def in nearsmth.c too) */
+
+
+#define XRES 100 /* Res of plots */
+
+/* Optional marker points for gamut mapping diagnosotic */
+struct {
+ int type; /* 1 = src point (xlate), 2 = dst point (no xlate) */
+ /* 0 = end marker */
+ double pos[3]; /* Position, (usually in Jab space) */
+ double col[3]; /* RGB color */
+} markers[] = {
+ { 0, }, /* End marker */
+ { 1, { 12.062, -0.87946, 0.97008 }, { 1.0, 0.3, 0.3 } }, /* Black point */
+ { 1, { 67.575411, -37.555250, -36.612862 }, { 1.0, 0.3, 0.3 } }, /* bad source in red (Red) */
+ { 1, { 61.003078, -44.466554, 1.922585 }, { 0.0, 1.0, 0.3 } }, /* good source in green */
+ { 2, { 49.294793, 50.749543, -51.383167 }, { 1.0, 0.0, 0.0 } },
+ { 2, { 42.783425, 49.089363, -37.823712 }, { 0.0, 1.0, 0.0 } },
+ { 2, { 41.222695, 63.911823, 37.695310 }, { 0.0, 1.0, 0.3 } }, /* destination in green */
+ { 1, { 41.951770, 60.220284, 34.788195 }, { 1.0, 0.3, 0.3 } }, /* source in red (Red) */
+ { 2, { 41.222695, 63.911823, 37.695310 }, { 0.3, 1.3, 0.3 } }, /* Dest in green */
+ { 1, { 85.117353, -60.807580, -22.195118 }, { 0.3, 0.3, 1 } }, /* Cyan Source (Blue) */
+ { 2, { 61.661622, -38.164411, -18.090824 }, { 1.0, 0.3, 0.3 } }, /* CMYK destination (Red) */
+ { 0 } /* End marker */
+};
+
+/* Optional marker rings for gamut mapping diagnosotic */
+struct {
+ int type; /* 1 = src ring point, 2 = ignore, */
+ /* 0 = end marker */
+ double ppoint[3]; /* Location of a point on the plane in source space */
+ double pnorm[3]; /* Plane normal direction in source space */
+ int nverts; /* Number of points to make ring */
+ double rad; /* Relative Radius from neutral to source surface (0.0 - 1.0) */
+ double scol[3]; /* Source RGB color */
+ double dcol[3]; /* Destination RGB color */
+} rings[] = {
+ { 0 }, /* End marker */
+ { 1,
+ { 60.0, 0.0, 0.0 }, { 1.0, 0.8, 0.0 }, /* plane point and normal */
+ 100, 1.0, /* 20 vertexes at source radius */
+ { 0.0, 1.0, 0.0 }, /* Green source */
+ { 1.0, 0.0, 0.0 } /* Red destination */
+ },
+ { 1,
+ { 60.0, 0.0, 0.0 }, { 1.0, 0.8, 0.0 }, /* plane point and normal */
+ 100, 0.9, /* 20 vertexes at source radius */
+ { 0.0, 1.0, 0.0 }, /* Green source */
+ { 1.0, 0.0, 0.0 } /* Red destination */
+ },
+ { 1,
+ { 60.0, 0.0, 0.0 }, { 1.0, 0.8, 0.0 }, /* plane point and normal */
+ 100, 0.8, /* 20 vertexes at source radius */
+ { 0.0, 1.0, 0.0 }, /* Green source */
+ { 1.0, 0.0, 0.0 } /* Red destination */
+ },
+ { 0 } /* End marker */
+};
+
+/* Degree to which the hue & saturation of the black point axes should be aligned: */
+#define GREYBPHSMF 0.0
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <fcntl.h>
+#include <string.h>
+#include <math.h>
+#include "counters.h"
+#include "icc.h"
+#include "numlib.h"
+#include "xicc.h"
+#include "gamut.h"
+#include "rspl.h"
+#include "gammap.h"
+#include "nearsmth.h"
+#include "vrml.h"
+#ifdef PLOT_LMAP
+#include "plot.h"
+#endif
+
+/* Callback context for enhancing the saturation of the clut values */
+typedef struct {
+ gamut *dst; /* Destination colorspace gamut */
+ double wp[3], bp[3];/* Destination colorspace white and black points */
+ double satenh; /* Saturation engancement value */
+} adjustsat;
+
+/* Callback context for making clut relative to white and black points */
+typedef struct {
+ double mat[3][4];
+} adjustwb;
+
+static void inv_grey_func(void *pp, double *out, double *in);
+static void adjust_wb_func(void *pp, double *out, double *in);
+static void adjust_sat_func(void *pp, double *out, double *in);
+
+#define XVRA 4.0 /* Extra mapping vertex ratio over no. tri verts from gamut */
+
+/* The smoothed near weighting control values. */
+/* These weightings setup the detailed behaviour of the */
+/* gamut mapping for the fully perceptual and saturation intents. */
+/* They are ordered here by increasing priority. A -ve value is ignored */
+
+/* Perceptual mapping weights, where smoothness and proportionality are important.. */
+gammapweights pweights[] = {
+ {
+ gmm_default, /* Non hue specific defaults */
+ { /* 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 */
+ },
+ 1.00 /* Chroma expansion 1 = none */
+ },
+ { /* 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.6, /* 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.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 */
+
+ 1.5, /* L error extra power with size, none = 1.0 */
+ 10.0 /* L error extra xover threshold in DE */
+ },
+ { /* Relative vector 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 ?) */
+ 10.0, /* Compression depth weight */
+ 10.0 /* Expansion depth weight */
+ },
+ {
+ 0.0 /* Fine tuning expansion weight, 0 - 1 */
+ }
+ },
+#ifdef NEVER
+ {
+ gmm_l_d_blue, /* Increase maintaining hue importance for blue */
+ {
+ {
+ -1.0, /* Cusp luminance alignment weighting 0 = none, 1 = full */
+ -1.0, /* Cusp chroma alignment weighting 0 = none, 1 = full */
+ 0.0 /* Cusp hue alignment weighting 0 = none, 1 = full */
+ },
+ -1.0 /* Chroma expansion 1 = none */
+ },
+ { /* Radial weighting */
+ -1.0, /* Radial error overall weight, 0 + */
+ -1.0, /* 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.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 */
+ -1.0, /* Black l dominance vs, c, 0 - 1 */
+
+ -1.0, /* White l threshold ratio to grey distance, 0 - 1 */
+ -1.0, /* Black l threshold ratio to grey distance, 0 - 1 */
+
+ -1.0, /* L error extra power, none = 1.0 */
+ -1.0 /* L error xover threshold in DE */
+ },
+ { /* Relative error preservation using smoothing */
+ -1.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 ?) */
+ -1.0, /* Compression depth weight */
+ -1.0 /* Expansion depth weight */
+ },
+ {
+ -1.0 /* Fine tuning expansion weight, 0 - 1 */
+ }
+ },
+#endif /* NEVER */
+ {
+ gmm_light_yellow, /* Treat yellow differently, to get purer result. */
+ {
+ {
+ 0.9, /* Cusp luminance alignment weighting 0 = none, 1 = full */
+ 0.8, /* Cusp chroma alignment weighting 0 = none, 1 = full */
+ 0.7 /* Cusp hue alignment weighting 0 = none, 1 = full */
+ },
+ 1.15 /* Chroma expansion 1 = none */
+ },
+ { /* Radial weighting */
+ -1.0, /* Radial error overall weight, 0 + */
+ -1.0, /* 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 */
+ -1.0, /* Hue dominance vs l+c, 0 - 1 */
+
+ -1.0, /* White l dominance vs, c, 0 - 1 */
+ -1.0, /* Grey l dominance vs, c, 0 - 1 */
+ -1.0, /* Black l dominance vs, c, 0 - 1 */
+
+ -1.0, /* White l threshold ratio to grey distance, 0 - 1 */
+ -1.0, /* Black l threshold ratio to grey distance, 0 - 1 */
+
+ -1.0, /* L error extra power, none = 1.0 */
+ -1.0 /* L error xover threshold in DE */
+ },
+ { /* Relative error preservation using smoothing */
+ 20.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 ?) */
+ -1.0, /* Compression depth weight */
+ -1.0 /* Expansion depth weight */
+ },
+ {
+ 0.5 /* Fine tuning expansion weight, 0 - 1 */
+ }
+ },
+ {
+ gmm_end,
+ }
+};
+double psmooth = 5.0; /* [5.0] Level of RSPL smoothing for perceptual, 1 = nominal */
+
+/* Saturation mapping weights, where saturation has priority over smoothness */
+gammapweights sweights[] = {
+ {
+ gmm_default, /* Non hue specific defaults */
+ { /* Cusp alignment control */
+ {
+ 0.6, /* Cusp luminance alignment weighting 0 = none, 1 = full */
+ 0.5, /* Cusp chroma alignment weighting 0 = none, 1 = full */
+ 0.6 /* Cusp hue alignment weighting 0 = none, 1 = full */
+ },
+ 1.05 /* Chroma expansion 1 = none */
+ },
+ { /* 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.4, /* Hue dominance vs l+c, 0 - 1 */
+
+ 0.6, /* White 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.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* */
+ },
+ { /* 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 ?) */
+ 10.0, /* Compression depth weight */
+ 10.0 /* Expansion depth weight */
+ },
+ {
+ 0.5 /* Fine tuning expansion weight, 0 - 1 */
+ }
+ },
+ {
+ gmm_light_yellow, /* Treat yellow differently, to get purer result. */
+ {
+ {
+ 1.0, /* Cusp luminance alignment weighting 0 = none, 1 = full */
+ 1.0, /* Cusp chroma alignment weighting 0 = none, 1 = full */
+ 1.0 /* Cusp hue alignment weighting 0 = none, 1 = full */
+ },
+ 1.20 /* Chroma expansion 1 = none */
+ },
+ { /* Radial weighting */
+ -1.0, /* Radial error overall weight, 0 + */
+ -1.0, /* 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.3, /* 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 */
+ -1.0, /* Black l dominance vs, c, 0 - 1 */
+
+ -1.0, /* White l threshold ratio to grey distance, 0 - 1 */
+ -1.0, /* Black l threshold ratio to grey distance, 0 - 1 */
+
+ -1.0, /* L error extra power, none = 1.0 */
+ -1.0 /* L error xover threshold in DE */
+ },
+ { /* Relative error preservation using 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. */
+ /* The depth is half the distance to the intersection of the */
+ /* vector to the other side of the gamut. (doesn't get triggered much ?) */
+ -1.0, /* Compression depth weight */
+ -1.0 /* Expansion depth weight */
+ },
+ {
+ -1.0 /* Fine tuning expansion weight, 0 - 1 */
+ }
+ },
+ {
+ gmm_end
+ }
+};
+double ssmooth = 2.0; /* Level of RSPL smoothing for saturation */
+
+/*
+ * Notes:
+ * The "knee" shape produced by the rspl (regular spline) code
+ * is not what one would expect for expansion. It is not
+ * symetrical with compression, and is less "sharp". This
+ * is due to the rspl "smoothness" criteria being based on
+ * grid value difference rather than smoothness being measured,
+ * as curvature. This means that the spline gets "stiffer" as
+ * it increases in slope.
+ * Possibly rspl could be improved in this respect ???
+ * (Doesn't matter for L compression now, because rspl is
+ * being inverted for expansion).
+ */
+
+static void del_gammap(gammap *s);
+static void domap(gammap *s, double *out, double *in);
+static void dopartialmap1(gammap *s, double *out, double *in);
+static void dopartialmap2(gammap *s, double *out, double *in);
+static gamut *parttransgamut(gammap *s, gamut *src);
+#ifdef PLOT_GAMUTS
+static void map_trans(void *cntx, double out[3], double in[3]);
+#endif
+
+/* Return a gammap to map from the input space to the output space */
+/* Return NULL on error. */
+gammap *new_gammap(
+ int verb, /* Verbose flag */
+ gamut *sc_gam, /* Source colorspace gamut */
+ 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 */
+ int dst_kbp, /* Use K only black point as dst gamut black point */
+ int dst_cmymap, /* masks C = 1, M = 2, Y = 4 to force 100% cusp map */
+ int rel_oride, /* 0 = normal, 1 = clip like, 2 = max relative */
+ int mapres, /* Gamut map resolution, typically 9 - 33 */
+ double *mn, /* If not NULL, set minimum mapping input range */
+ double *mx, /* for rspl grid. */
+ char *diagname /* If non-NULL, write a gamut mapping diagnostic WRL */
+) {
+ gmm_BPmap bph = gmm_bendBP; /* Prefered algorithm */
+// gmm_BPmap bph = gmm_clipBP; /* Alternatives tried */
+// gmm_BPmap bph = gmm_BPadpt;
+// gmm_BPmap bph = gmm_noBPadpt;
+
+ 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 */
+ double s_ga_wp[3]; /* Source (image) gamut white point */
+ double s_ga_bp[3]; /* Source (image) gamut black point */
+ double d_cs_wp[3]; /* Destination colorspace white point */
+ double d_cs_bp[3]; /* Destination colorspace black point */
+
+ double sr_cs_wp[3]; /* Source rotated colorspace white point */
+ double sr_cs_bp[3]; /* Source rotated colorspace black point */
+ double sr_ga_wp[3]; /* Source rotated (image) gamut white point */
+ double sr_ga_bp[3]; /* Source rotated (image) gamut black point */
+ double dr_cs_wp[3]; /* Target (gmi->greymf aligned) white point */
+ double dr_cs_bp[3]; /* Target (gmi->greymf aligned) black point */
+ double dr_be_bp[3]; /* Bend at start in source neutral axis direction */
+ /* Target black point (Same as dr_cs_bp[] otherwise) */
+
+ double sl_cs_wp[3]; /* Source rotated and L mapped colorspace white point */
+ double sl_cs_bp[3]; /* Source rotated and L mapped colorspace black point */
+
+ double s_mt_wp[3]; /* Overall source mapping target white point (used for finetune) */
+ double s_mt_bp[3]; /* Overall source mapping target black point (used for finetune) */
+ 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) */
+
+ 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 */
+ int ngreyp = 0; /* Number of grey axis mapping points */
+ int ngamp = 0; /* Number of gamut mapping points */
+ double xvra = XVRA; /* Extra ss vertex ratio to src gamut vertex count */
+ int j;
+
+#if defined(PLOT_LMAP) || defined(PLOT_GAMUTS) || defined(PLOT_3DKNEES)
+ fprintf(stderr,"##### A gammap.c PLOT is #defined ####\n");
+#endif
+
+ if (verb) {
+ xicc_dump_gmi(gmi);
+ printf("Gamut map resolution: %d\n",mapres);
+ if (si_gam != NULL)
+ printf("Image gamut supplied\n");
+ switch(bph) {
+ case gmm_clipBP: printf("Neutral axis no-adapt extend and clip\n"); break;
+ case gmm_BPadpt: printf("Neutral axis fully adapt\n"); break;
+ case gmm_bendBP: printf("Neutral axis no-adapt extend and bend\n"); break;
+ case gmm_noBPadpt: printf("Neutral axis no-adapt\n"); break;
+ }
+ }
+
+ /* Allocate the object */
+ if ((s = (gammap *)calloc(1, sizeof(gammap))) == NULL)
+ error("gammap: calloc failed on gammap object");
+
+ /* Setup methods */
+ s->del = del_gammap;
+ s->domap = domap;
+
+ /* Now create everything */
+
+ /* Grab the white and black points */
+ if (src_kbp) {
+ // ~~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, 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;
+ }
+ }
+
+ /* If source space is source gamut */
+ 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 sourcegamut */
+ } else {
+
+ 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");
+ free(s);
+ return NULL;
+ }
+ } else {
+ if (si_gam->getwb(si_gam, NULL, NULL, NULL, s_ga_wp, s_ga_bp, NULL)) {
+ fprintf(stderr,"gamut map: Unable to read source gamut white and black points\n");
+ free(s);
+ 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");
+ 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");
+ free(s);
+ return NULL;
+ }
+ }
+
+#ifdef VERBOSE
+ if (verb) {
+ if (src_kbp)
+ printf("Using Src K only black point\n");
+
+ if (dst_kbp)
+ printf("Using Dst K only black point\n");
+
+ printf("Src colorspace white/black are %f %f %f, %f %f %f\n",
+ s_cs_wp[0], s_cs_wp[1], s_cs_wp[2], s_cs_bp[0], s_cs_bp[1], s_cs_bp[2]);
+
+ printf("Src gamut white/black are %f %f %f, %f %f %f\n",
+ s_ga_wp[0], s_ga_wp[1], s_ga_wp[2], s_ga_bp[0], s_ga_bp[1], s_ga_bp[2]);
+
+ printf("Dst colorspace white/black are %f %f %f, %f %f %f\n",
+ d_cs_wp[0], d_cs_wp[1], d_cs_wp[2], d_cs_bp[0], d_cs_bp[1], d_cs_bp[2]);
+ }
+#endif /* VERBOSE */
+
+ /* ------------------------------------ */
+ /* Figure out the destination grey axis alignment */
+ /* This is all done using colorspace white & black points */
+ {
+ double t, svl, dvl;
+ double wrot[3][3]; /* Rotation about 0,0,0 to match white points */
+ double sswp[3], ssbp[3]; /* Temporary source white & black points */
+ double fawp[3], fabp[3]; /* Fully adapted destination white & black */
+ double hawp[3], habp[3]; /* Half (full white, not black) adapted destination w & b */
+
+ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+ /* The first task is to decide what our target destination */
+ /* white and black points are going to be. */
+
+ /* Figure out what our initial target destination white point is going to be: */
+
+ /* Compute source white and black points with same L value as the destination */
+ t = (d_cs_wp[0] - s_cs_bp[0])/(s_cs_wp[0] - s_cs_bp[0]);
+ for (j = 0; j < 3; j++)
+ sswp[j] = s_cs_bp[j] + t * (s_cs_wp[j] - s_cs_bp[j]);
+
+ t = (d_cs_bp[0] - s_cs_wp[0])/(s_cs_bp[0] - s_cs_wp[0]);
+ for (j = 0; j < 3; j++)
+ ssbp[j] = s_cs_wp[j] + t * (s_cs_bp[j] - s_cs_wp[j]);
+
+ /* The raw grey axis alignment target is a blend between the */
+ /* source colorspace (NOT gamut) and the destination */
+ /* colorspace. */
+
+ for (j = 0; j < 3; j++) {
+ dr_cs_wp[j] = gmi->greymf * d_cs_wp[j] + (1.0 - gmi->greymf) * sswp[j];
+ dr_cs_bp[j] = gmi->greymf * d_cs_bp[j] + (1.0 - gmi->greymf) * ssbp[j];
+ }
+
+#ifdef VERBOSE
+ if (verb) {
+ printf("Target (blended) dst wp/bp = %f %f %f, %f %f %f\n",
+ dr_cs_wp[0], dr_cs_wp[1], dr_cs_wp[2], dr_cs_bp[0], dr_cs_bp[1], dr_cs_bp[2]);
+ }
+#endif /* VERBOSE */
+
+ /* Compute full adaptation target destinations */
+ for (j = 0; j < 3; j++) {
+ fawp[j] = dr_cs_wp[j]; /* White fully adapted */
+ fabp[j] = dr_cs_bp[j]; /* Black fully adapted */
+ }
+
+ /* Clip the target grey axis to the destination gamut */
+ if (d_gam->vector_isect(d_gam, fabp, fawp, fabp, fawp, NULL, NULL, NULL, NULL) == 0)
+ error("gamut: vector_isect failed!");
+
+ /* To work around the problem that vector_isect() is not entirely accurate, */
+ /* special case the situation where gmi->greymf == 1.0 */
+ if (gmi->greymf > 0.99) {
+ for (j = 0; j < 3; j++) {
+ fawp[j] = d_cs_wp[j];
+ fabp[j] = d_cs_bp[j];
+ }
+ }
+
+ /* If dst_kbp is set, then clipping to the dest gamut doesn't do what we want, */
+ /* since it extends the black to a full composite black point. */
+ /* A "K only" gamut is hard to define, so do a hack: */
+ /* scale fabp[] towards fawp[] so that it has the same L as */
+ /* the destination K only black point. */
+ if (dst_kbp && fabp[0] < d_cs_bp[0]) {
+ t = (d_cs_bp[0] - fawp[0])/(fabp[0] - fawp[0]);
+
+ for (j = 0; j < 3; j++)
+ fabp[j] = fawp[j] + t * (fabp[j] - fawp[j]);
+ }
+
+ /* Compute half adapted (full white, not black) target destinations */
+ for (j = 0; j < 3; j++)
+ hawp[j] = dr_cs_wp[j]; /* White fully adapted */
+
+ /* Compute the rotation matrix that maps the source white point */
+ /* onto the target white point. */
+ icmRotMat(wrot, sswp, dr_cs_wp);
+
+ /* Compute the target black point as the rotated source black point */
+ icmMulBy3x3(habp, wrot, s_cs_bp);
+
+ /* Now intersect the target white and black points with the destination */
+ /* colorspace gamut to arrive at the best possible in gamut values for */
+ /* the target white and black points. */
+ if (d_gam->vector_isect(d_gam, habp, hawp, habp, hawp, NULL, NULL, NULL, NULL) == 0)
+ error("gamut: vector_isect failed!");
+
+ /* To work around the problem that vector_isect() is not entirely accurate, */
+ /* special case the situation where gmi->greymf == 1.0 */
+ if (gmi->greymf > 0.99) {
+ for (j = 0; j < 3; j++) {
+ hawp[j] = d_cs_wp[j];
+ }
+ }
+
+ /* If dst_kbp is set, then clipping to the dest gamut doesn't do what we want, */
+ /* since it extends the black to a full composite black point. */
+ /* A "K only" gamut is hard to define, so do a hack: */
+ /* scale habp[] towards hawp[] so that it has the same L as */
+ /* the destination K only black point. */
+ if (dst_kbp && habp[0] < d_cs_bp[0]) {
+ t = (d_cs_bp[0] - hawp[0])/(habp[0] - hawp[0]);
+
+ for (j = 0; j < 3; j++)
+ habp[j] = hawp[j] + t * (habp[j] - hawp[j]);
+ }
+
+ /* Now decide the detail of the white and black alignment */
+ if (bph == gmm_BPadpt || bph == gmm_bendBP) { /* Adapt to destination white and black */
+
+
+ /* Use the fully adapted white and black points */
+ for (j = 0; j < 3; j++) {
+ dr_cs_wp[j] = fawp[j];
+ dr_cs_bp[j] = fabp[j];
+ }
+
+ if (bph == gmm_bendBP) {
+
+ /* Extend the half adapted (white = dst, black = src) black point */
+ /* to the same L as the target (dst), to use as the initial (bent) black point */
+ t = (dr_cs_bp[0] - dr_cs_wp[0])/(habp[0] - dr_cs_wp[0]);
+ for (j = 0; j < 3; j++)
+ dr_be_bp[j] = dr_cs_wp[j] + t * (habp[j] - dr_cs_wp[j]);
+
+ } else {
+
+ /* Set bent black point target to be the same as our actual */
+ /* black point target, so that the "bend" code does nothing. */
+ for (j = 0; j < 3; j++)
+ dr_be_bp[j] = dr_cs_bp[j];
+ }
+
+ } else { /* Adapt to destination white but not black */
+
+ /* Use the half adapted (white = dst, black = src) white and black points */
+ for (j = 0; j < 3; j++) {
+ dr_cs_wp[j] = hawp[j];
+ dr_cs_bp[j] = habp[j];
+ }
+
+#ifdef VERBOSE
+ if (verb) {
+ printf("Adapted target wp/bp = %f %f %f, %f %f %f\n",
+ dr_cs_wp[0], dr_cs_wp[1], dr_cs_wp[2], dr_cs_bp[0], dr_cs_bp[1], dr_cs_bp[2]);
+ }
+#endif
+ if (bph == gmm_clipBP) {
+
+ /* Extend the target black point to accomodate the */
+ /* bent or clipped destination space L* range */
+ if (fabp[0] < dr_cs_bp[0]) {
+ t = (fabp[0] - dr_cs_wp[0])/(dr_cs_bp[0] - dr_cs_wp[0]);
+ for (j = 0; j < 3; j++)
+ dr_cs_bp[j] = dr_cs_wp[j] + t * (dr_cs_bp[j] - d_cs_wp[j]);
+ }
+ }
+
+ /* Set the bent black point target to be the same as our actual */
+ /* black point target, so that the "bend" code does nothing. */
+ for (j = 0; j < 3; j++)
+ dr_be_bp[j] = dr_cs_bp[j];
+ }
+
+#ifdef VERBOSE
+ if (verb) {
+ printf("Adapted & extended tgt wp/bp = %f %f %f, %f %f %f\n",
+ dr_cs_wp[0], dr_cs_wp[1], dr_cs_wp[2], dr_cs_bp[0], dr_cs_bp[1], dr_cs_bp[2]);
+ }
+#endif /* VERBOSE */
+
+ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+ /* Now we need to figure out what origin alignment is needed, as well as */
+ /* making sure the vectors are the same length to avoid rescaling. */
+ /* (Scaling is meant to be done with the L curve though.) */
+
+ /* Create temporary source white point that has the same L as the */
+ /* target destination white point. */
+ t = (dr_cs_wp[0] - s_cs_bp[0])/(s_cs_wp[0] - s_cs_bp[0]);
+ for (j = 0; j < 3; j++)
+ sswp[j] = s_cs_bp[j] + t * (s_cs_wp[j] - s_cs_bp[j]);
+
+ /* Create temporary source black point that will form a vector to the src white */
+ /* point with same length as the target destination black->white vector. */
+ for (svl = dvl = 0.0, j = 0; j < 3; j++) {
+ double tt;
+ tt = sswp[j] - s_cs_bp[j];
+ svl += tt * tt;
+ tt = dr_cs_wp[j] - dr_cs_bp[j];
+ dvl += tt * tt;
+ }
+ svl = sqrt(svl);
+ dvl = sqrt(dvl);
+ for (j = 0; j < 3; j++)
+ ssbp[j] = sswp[j] + dvl/svl * (s_cs_bp[j] - sswp[j]);
+
+#ifdef VERBOSE
+ if (verb) {
+ printf("Rotate matrix src wp/bp = %f %f %f, %f %f %f\n",
+ sswp[0], sswp[1], sswp[2], ssbp[0], ssbp[1], ssbp[2]);
+ printf("Rotate matrix dst wp/bp = %f %f %f, %f %f %f\n",
+ dr_cs_wp[0], dr_cs_wp[1], dr_cs_wp[2], dr_cs_bp[0], dr_cs_bp[1], dr_cs_bp[2]);
+ }
+#endif /* VERBOSE */
+
+ /* Now create the general rotation and translation to map the source grey */
+ /* axis to our destination grey axis. */
+ icmVecRotMat(s->grot, sswp, ssbp, dr_cs_wp, dr_cs_bp);
+
+ /* And create the inverse as well: */
+ icmVecRotMat(s->igrot, dr_cs_wp, dr_cs_bp, sswp, ssbp);
+
+ /* Create rotated versions of source colorspace & image white and */
+ /* black points for use from now on, given that rotation will */
+ /* be applied first to all source points. */
+ icmMul3By3x4(sr_cs_wp, s->grot, s_cs_wp);
+ icmMul3By3x4(sr_cs_bp, s->grot, s_cs_bp);
+ icmMul3By3x4(sr_ga_wp, s->grot, s_ga_wp);
+ icmMul3By3x4(sr_ga_bp, s->grot, s_ga_bp);
+
+#ifdef VERBOSE
+ if (verb) {
+ printf("Bend target bp = %f %f %f\n",
+ dr_be_bp[0], dr_be_bp[1], dr_be_bp[2]);
+ printf("Rotated source grey axis wp/bp %f %f %f, %f %f %f\n",
+ sr_cs_wp[0], sr_cs_wp[1], sr_cs_wp[2], sr_cs_bp[0], sr_cs_bp[1], sr_cs_bp[2]);
+ printf("Rotated gamut grey axis wp/bp %f %f %f, %f %f %f\n",
+ sr_ga_wp[0], sr_ga_wp[1], sr_ga_wp[2], sr_ga_bp[0], sr_ga_bp[1], sr_ga_bp[2]);
+ printf("Destination axis target wp/bp %f %f %f, %f %f %f\n",
+ dr_cs_wp[0], dr_cs_wp[1], dr_cs_wp[2], dr_cs_bp[0], dr_cs_bp[1], dr_cs_bp[2]);
+ }
+#endif
+ }
+
+#ifdef NEVER
+sr_cs_wp[0] = 100.0;
+sr_cs_bp[0] = 30.0;
+dr_cs_wp[0] = 80.0;
+dr_cs_bp[0] = 10.0;
+glumknf = 1.0;
+#endif /* NEVER */
+
+ /* Create the mapping points needed to build the 1D L mapping rspl. */
+ /* 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 a gamut/image range. */
+ {
+ double swL, dwL; /* Source and destination white point L */
+ double sbL, dbL; /* Source and destination black point L */
+ int j;
+ double t;
+
+ /* Setup white point mapping */
+ if (sr_cs_wp[0] <= dr_cs_wp[0]) { /* Needs possible expansion */
+ swL = sr_cs_wp[0];
+ dwL = gmi->glumwexf * dr_cs_wp[0] + (1.0 - gmi->glumwexf) * sr_cs_wp[0];
+
+ } else {
+ 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_ga_wp[0];
+ dwL = dr_cs_wp[0];
+
+ } else { /* Neither needed */
+ swL = sr_ga_wp[0];
+ dwL = sr_ga_wp[0];
+ }
+ }
+
+ /* Setup black point mapping */
+ if (sr_cs_bp[0] >= dr_cs_bp[0]) { /* Needs possible expansion */
+ sbL = sr_cs_bp[0];
+ dbL = gmi->glumbexf * dr_cs_bp[0] + (1.0 - gmi->glumbexf) * sr_cs_bp[0];
+
+ } else {
+ 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_ga_bp[0];
+ dbL = dr_cs_bp[0];
+
+ } else { /* Neither needed */
+ sbL = sr_ga_bp[0];
+ dbL = sr_ga_bp[0];
+ }
+ }
+
+ /* To ensure symetry between compression and expansion, always create RSPL */
+ /* for compression and its inverse, and then swap grey and igrey rspl to compensate. */
+ if ((dwL - dbL) > (swL - sbL))
+ revrspl = 1;
+
+ /* White point end */
+ lpnts[ngreyp].p[0] = swL;
+ lpnts[ngreyp].v[0] = dwL;
+ lpnts[ngreyp++].w = 10.0; /* Must go through here */
+
+ /* Black point end */
+ lpnts[ngreyp].p[0] = sbL;
+ lpnts[ngreyp].v[0] = dbL;
+ lpnts[ngreyp++].w = 10.0; /* Must go through here */
+
+//printf("~1 white loc %f, val %f\n",swL,dwL);
+//printf("~1 black loc %f, val %f\n",sbL,dbL);
+
+#ifdef USE_GLUMKNF
+ if (gmi->glumknf < 0.05)
+#endif /* USE_GLUMKNF */
+ { /* 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 cplv; /* Center point location and value */
+ double kppos = 0.30; /* Knee point ratio between white/black & center */
+ double kwl, kbl, kwv, kbv; /* Knee point values and locations */
+ double kwx, kbx; /* Knee point extra */
+
+
+//printf("sbL = %f, swL = %f\n",sbL,swL);
+//printf("dbL = %f, dwL = %f\n",dbL,dwL);
+
+ /* Center point */
+ cplv = cppos * (swL - sbL) + sbL;
+//printf("~1 computed cplv = %f\n",cplv);
+
+#ifdef NEVER /* Don't use a center point */
+ lpnts[ngreyp].p[0] = cplv;
+ lpnts[ngreyp].v[0] = cplv;
+ lpnts[ngreyp++].w = 0.5;
+#endif
+
+//printf("~1 black half diff = %f\n",dbL - sbL);
+//printf("~1 white half diff = %f\n",dwL - swL);
+
+ /* Knee point locations */
+ kwl = kppos * (cplv - swL) + swL;
+ kbl = kppos * (cplv - sbL) + sbL;
+
+ /* Extra compression for white and black knees */
+ 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;
+
+ /* Knee point values */
+ kwv = (dwL + kwx - cplv) * (kwl - cplv)/(swL - cplv) + cplv;
+ if (kwv > dwL) /* Sanity check */
+ kwv = dwL;
+
+ kbv = (dbL - kbx - cplv) * (kbl - cplv)/(sbL - cplv) + cplv;
+ if (kbv < dbL) /* Sanity check */
+ kbv = dbL;
+
+
+//printf("~1 kbl = %f, kbv = %f\n",kbl, kbv);
+//printf("~1 kwl = %f, kwv = %f\n",kwl, kwv);
+
+ /* Emphasise points to cause "knee" curve */
+ lpnts[ngreyp].p[0] = kwl;
+ lpnts[ngreyp].v[0] = kwv;
+ lpnts[ngreyp++].w = gmi->glumknf * gmi->glumknf;
+
+ lpnts[ngreyp].p[0] = kbl;
+ lpnts[ngreyp].v[0] = kbv;
+ lpnts[ngreyp++].w = 1.5 * gmi->glumknf * 1.5 * gmi->glumknf;
+ }
+#endif /* USE_GLUMKNF */
+
+ /* Remember our source and destinatio mapping targets */
+ /* so that we can use them for fine tuning later. */
+
+ /* We scale the source and target white and black */
+ /* points to match the L values of the source and destination */
+ /* L curve mapping, as this is how we have chosen the */
+ /* white and black point mapping for the link. */
+ /* Put them back in pre-rotated space, so that we can */
+ /* check the overall transform of the white and black points. */
+ t = (swL - sr_cs_bp[0])/(sr_cs_wp[0] - sr_cs_bp[0]);
+ for (j = 0; j < 3; j++)
+ s_mt_wp[j] = sr_cs_bp[j] + t * (sr_cs_wp[j] - sr_cs_bp[j]);
+ icmMul3By3x4(s_mt_wp, s->igrot, s_mt_wp);
+
+ t = (sbL - sr_cs_wp[0])/(sr_cs_bp[0] - sr_cs_wp[0]);
+ for (j = 0; j < 3; j++)
+ s_mt_bp[j] = sr_cs_wp[j] + t * (sr_cs_bp[j] - sr_cs_wp[j]);
+//printf("~1 check black point rotated = %f %f %f\n",s_mt_bp[0],s_mt_bp[1],s_mt_bp[2]);
+ icmMul3By3x4(s_mt_bp, s->igrot, s_mt_bp);
+//printf("~1 check black point prerotated = %f %f %f\n",s_mt_bp[0],s_mt_bp[1],s_mt_bp[2]);
+
+ t = (dwL - dr_cs_bp[0])/(dr_cs_wp[0] - dr_cs_bp[0]);
+ for (j = 0; j < 3; j++)
+ d_mt_wp[j] = dr_cs_bp[j] + t * (dr_cs_wp[j] - dr_cs_bp[j]);
+
+ for (j = 0; j < 3; j++)
+ d_mt_bp[j] = dr_cs_wp[j] + t * (dr_cs_bp[j] - dr_cs_wp[j]);
+ }
+
+ /* We now create the 1D rspl L map, that compresses or expands the luminence */
+ /* 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;
+
+ if (revrspl) { /* Invert creation and usage for symetry between compress and exp. */
+ int i;
+ for (i = 0; i < ngreyp; i++) {
+ double tt = lpnts[i].p[0]; /* Swap source and dest */
+ lpnts[i].p[0] = lpnts[i].v[0];
+ lpnts[i].v[0] = tt;
+ }
+ }
+
+ /* Create a 1D rspl, that is used to */
+ /* form the overall L compression mapping. */
+ if ((s->grey = new_rspl(RSPL_NOFLAGS, 1, 1)) == NULL) /* Allocate 1D -> 1D */
+ 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);
+ }
+#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");
+ }
+
+ /* 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 (revrspl) { /* Swap to compensate for expansion */
+ rspl *tt = s->grey;
+ s->grey = s->igrey;
+ s->igrey = tt;
+ }
+ }
+
+#ifdef PLOT_LMAP
+ { /* Plot the 1D mapping */
+ double xx[XRES];
+ double y1[XRES];
+ int i;
+
+ for (i = 0; i < XRES; i++) {
+ double x;
+ co cp; /* Conversion point */
+ x = sr_cs_bp[0] + (i/(double)(XRES-1)) * (sr_cs_wp[0] - sr_cs_bp[0]);
+ xx[i] = x;
+ cp.p[0] = x;
+ s->grey->interp(s->grey, &cp);
+ y1[i] = cp.v[0];
+ }
+ do_plot(xx,y1,NULL,NULL,XRES);
+ }
+#endif /* PLOT_LMAP */
+
+ {
+ /* We want to rotate and then map L independently of everything else, */
+ /* so transform source csape & image gamuts through the rotation and L mapping */
+ /* before we create the surface 3D mapping from them */
+
+ /* Create L mapped versions of rotated src colorspace white/black points */
+#ifdef NEVER
+ co cp;
+ double t;
+ int i;
+
+ cp.p[0] = sr_cs_wp[0];
+ s->grey->interp(s->grey, &cp);
+
+ t = (cp.v[0] - sr_cs_bp[0])/(sr_cs_wp[0] - sr_cs_bp[0]);
+ for (j = 0; j < 3; j++)
+ sl_cs_wp[j] = sr_cs_bp[j] + t * (sr_cs_wp[j] - sr_cs_bp[j]);
+
+ cp.p[0] = sr_cs_bp[0];
+ s->grey->interp(s->grey, &cp);
+ t = (cp.v[0] - sr_cs_wp[0])/(sr_cs_bp[0] - sr_cs_wp[0]);
+ for (j = 0; j < 3; j++)
+ sl_cs_bp[j] = sr_cs_wp[j] + t * (sr_cs_bp[j] - sr_cs_wp[j]);
+#else
+ dopartialmap1(s, sl_cs_wp, s_cs_wp);
+ dopartialmap1(s, sl_cs_bp, s_cs_bp);
+#endif
+
+#ifdef VERBOSE
+ if (verb) {
+ printf("Mapped source grey axis wp/bp %f %f %f, %f %f %f\n",
+ sl_cs_wp[0], sl_cs_wp[1], sl_cs_wp[2], sl_cs_bp[0], sl_cs_bp[1], sl_cs_bp[2]);
+ }
+#endif
+
+ if ((scl_gam = parttransgamut(s, sc_gam)) == NULL) {
+ fprintf(stderr,"gamut map: parttransgamut failed\n");
+ free(s);
+ return NULL;
+ }
+
+ if (sc_gam == si_gam)
+ sil_gam = scl_gam;
+
+ else {
+ if ((sil_gam = parttransgamut(s, si_gam)) == NULL) {
+ fprintf(stderr,"gamut map: parttransgamut failed\n");
+ free(s);
+ return NULL;
+ }
+ }
+ }
+
+ /* Create all the 3D->3D gamut mapping points and 3D rspl, */
+ /* if there is any compression or expansion to do. */
+ if (gmi->gamcpf > 1e-6 || gmi->gamexf > 1e-6) {
+ cow *gpnts = NULL; /* Mapping points to create gamut mapping */
+ int nspts; /* Number of source gamut surface points */
+ int rgridpts; /* Number of range surface grid points */
+ int i, j;
+ datai il, ih;
+ datao ol, oh;
+ int gres[MXDI];
+ double avgdev[MXDO];
+ 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], 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) */
+ cgats *locus = NULL; /* Diagnostic locus to plot in wrl, NULL if none */
+
+#ifdef PLOT_3DKNEES
+typedef struct {
+ double v0[3], v1[3];
+} p3dk_lpoint;
+ p3dk_lpoint *p3dk_locus;
+ int p3dk_ix = 0;
+#endif /* PLOT_3DKNEES */
+
+ /* Get the maximum number of points that will be created */
+ nspts = near_smooth_np(scl_gam, sil_gam, d_gam, xvra);
+
+ rgridpts = 0;
+#ifdef USE_BOUND
+ if (defrgrid >= 2) {
+ rgridpts = defrgrid * defrgrid * defrgrid
+ - (defrgrid -2) * (defrgrid -2) * (defrgrid -2);
+ }
+#endif
+
+ if ((gpnts = (cow *)malloc((nres + 3 * nspts + rgridpts) * sizeof(cow))) == NULL) {
+ fprintf(stderr,"gamut map: Malloc of mapping setup points failed\n");
+ s->grey->del(s->grey);
+ s->igrey->del(s->igrey);
+ if (sil_gam != scl_gam)
+ sil_gam->del(sil_gam);
+ scl_gam->del(scl_gam);
+ free(s);
+ return NULL;
+ }
+
+#ifdef PLOT_3DKNEES
+ if ((p3dk_locus = (p3dk_lpoint *)malloc((2 * nspts) * sizeof(p3dk_lpoint))) == NULL)
+ error("gamut: Diagnostic array p3dk_locus malloc failed");
+#endif /* PLOT_3DKNEES */
+
+ /* ------------------------------------------- */
+ /* Finish off the grey axis mapping by creating the */
+ /* grey axis 3D->3D mapping points */
+ /* We use 4 times the grid density, and create */
+ /* points that span the source colorspace (this may exceed) */
+ /* the source image gamut, and map to points outside the */
+ /* destination gamut) */
+
+ /* See how much to bend the black - compute the color difference */
+ /* We start out in the direction of dr_be_bp at white, and at */
+ /* the end we bend towards the overall bp dr_cs_bp */
+ /* (brad will be 0 for non gmm_bendBP because dr_be_bp dr_cs_bp */
+ for (brad = 0.0, i = 1; i < 3; i++) {
+ double tt = dr_be_bp[i] - dr_cs_bp[i];
+ brad += tt * tt;
+ }
+ brad = sqrt(brad);
+
+//printf("~1 brad = %f, Bend target = %f %f %f, straight = %f %f %f\n",
+//brad, dr_be_bp[0], dr_be_bp[1], dr_be_bp[2], dr_cs_bp[0], dr_cs_bp[1], dr_cs_bp[2]);
+
+#ifdef USE_GREYMAP
+ for (i = 0; i < nres; i++) { /* From black to white */
+ double t;
+ double bv[3]; /* Bent (initial) destination value */
+ double dv[3]; /* Straight (final) destination value */
+ double wt = 1.0; /* Default grey axis point weighting */
+
+ /* Create source grey axis point */
+ t = i/(nres - 1.0);
+
+ /* Cover L = 0.0 to 100.0 */
+ t = ((100.0 * t) - sl_cs_bp[0])/(sl_cs_wp[0] - sl_cs_bp[0]);
+ for (j = 0; j < 3; j++)
+ gpnts[ngamp].p[j] = sl_cs_bp[j] + t * (sl_cs_wp[j] - sl_cs_bp[j]);
+
+ /* L values are the same, as they have been mapped prior to 3D */
+ gpnts[ngamp].v[0] = gpnts[ngamp].p[0];
+
+ /* Figure destination point on initial bent grey axis */
+ t = (gpnts[ngamp].v[0] - dr_cs_wp[0])/(dr_be_bp[0] - dr_cs_wp[0]);
+ for (j = 0; j < 3; j++)
+ bv[j] = dr_cs_wp[j] + t * (dr_be_bp[j] - dr_cs_wp[j]);
+//printf("~1 t = %f, bent dest %f %f %f\n",t, bv[0], bv[1],bv[2]);
+
+ /* Figure destination point on final straight grey axis */
+ t = (gpnts[ngamp].v[0] - dr_cs_wp[0])/(dr_cs_bp[0] - dr_cs_wp[0]);
+ for (j = 0; j < 3; j++)
+ dv[j] = dr_cs_wp[j] + t * (dr_cs_bp[j] - dr_cs_wp[j]);
+//printf("~1 t = %f, straight dest %f %f %f\n",t, dv[0], dv[1],dv[2]);
+
+ /* Figure out a blend value between the bent value */
+ /* and the straight value, so that it curves smoothly from */
+ /* one to the other. */
+ if (brad > 0.001) {
+ double ty;
+ t = ((dr_cs_bp[0] + brad) - gpnts[ngamp].v[0])/brad;
+ if (t < 0.0)
+ t = 0.0;
+ else if (t > 1.0)
+ t = 1.0;
+ /* Make it a spline ? */
+ t = t * t * (3.0 - 2.0 * t);
+ ty = t * t * (3.0 - 2.0 * t); /* spline blend value */
+ t = (1.0 - t) * ty + t * t; /* spline at t == 0, linear at t == 1 */
+
+ wt *= (1.0 + t * brad); /* Increase weigting with the bend */
+
+ } else {
+ t = 0.0; /* stick to straight, it will be close anyway. */
+ }
+
+ for (j = 0; j < 3; j++) /* full straight when t == 1 */
+ gpnts[ngamp].v[j] = t * dv[j] + (1.0 - t) * bv[j];
+ 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 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],
+ gpnts[ngamp].w);
+#endif
+ ngamp++;
+ }
+#endif /* USE_GREYMAP */
+
+ /* ---------------------------------------------------- */
+ /* Do preliminary computation of the rspl input and output bounding values */
+ for (j = 0; j < 3; j++) {
+ il[j] = ol[j] = 1e60;
+ ih[j] = oh[j] = -1e60;
+ }
+
+ /* From grey axis points */
+ for (i = 0; i < ngamp; i++) {
+ for (j = 0; j < 3; j++) {
+ if (gpnts[i].p[j] < il[j])
+ il[j] = gpnts[i].p[j];
+ if (gpnts[i].p[j] > ih[j])
+ ih[j] = gpnts[i].p[j];
+ }
+ }
+
+ /* From the source gamut */
+ {
+ double tmx[3], tmn[3];
+ scl_gam->getrange(scl_gam, tmn, tmx);
+ for (j = 0; j < 3; j++) {
+ if (tmn[j] < il[j])
+ il[j] = tmn[j];
+ if (tmx[j] > ih[j])
+ ih[j] = tmx[j];
+ }
+ }
+
+ /* from input arguments override */
+ if (mn != NULL && mx != NULL) {
+
+ for (j = 0; j < 3; j++) {
+ if (mn[j] < il[j])
+ il[j] = mn[j];
+ if (mx[j] > ih[j])
+ ih[j] = mx[j];
+ }
+ }
+
+ /* From the destination gamut */
+ {
+ double tmx[3], tmn[3];
+ d_gam->getrange(d_gam, tmn, tmx);
+ for (j = 0; j < 3; j++) {
+ if (tmn[j] < ol[j])
+ ol[j] = tmn[j];
+ if (tmx[j] > oh[j])
+ oh[j] = tmx[j];
+ }
+ }
+
+ /* ---------------------------------------------------- */
+ /* Deal with gamut hull guide vector creation. */
+
+ /* For compression, create a mapping for each vertex of */
+ /* the source gamut (image) surface towards the destination gamut */
+ /* For expansion, do the opposite. */
+
+ /* Convert from compact to explicit hextant weightings */
+ if (expand_weights(xpweights, pweights)
+ || expand_weights(xsweights, sweights)) {
+ fprintf(stderr,"gamut map: expand_weights() failed\n");
+ s->grey->del(s->grey);
+ s->igrey->del(s->igrey);
+ if (sil_gam != scl_gam)
+ sil_gam->del(sil_gam);
+ scl_gam->del(scl_gam);
+ free(s);
+ return NULL;
+ }
+ /* 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) {
+ tweak_weights(xwh, dst_cmymap, rel_oride);
+ }
+
+ /* Create the near point mapping, which is our fundamental gamut */
+ /* hull to gamut hull mapping. */
+ 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, il, ih, ol, oh);
+ if (nsm == NULL) {
+ fprintf(stderr,"Creating smoothed near points failed\n");
+ s->grey->del(s->grey);
+ s->igrey->del(s->igrey);
+ if (sil_gam != scl_gam)
+ sil_gam->del(sil_gam);
+ scl_gam->del(scl_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 */
+ /* as we create the final 3D gamut mapping rspl */
+ /* (The plot is of the already rotated and L mapped source space) */
+ {
+ int doaxes = 0;
+
+#ifdef PLOT_AXES
+ doaxes = 1;
+#endif
+ if (diagname != NULL)
+ wrl = new_vrml(diagname, doaxes, 0);
+#ifdef PLOT_DIAG_WRL
+ else
+ wrl = new_vrml("gammap.wrl", doaxes, 0);
+#endif
+ }
+
+ if (wrl != NULL) {
+ /* See if there is a diagnostic locus to plot too */
+ if ((locus = new_cgats()) == NULL)
+ error("Failed to create cgats object");
+
+ locus->add_other(locus, "TS");
+
+ if (locus->read_name(locus, "locus.ts")) {
+ locus->del(locus);
+ locus = NULL;
+ } else {
+ if (verb)
+ printf("!! Found diagnostic locus.ts file !!\n");
+ /* locus will be added later */
+ }
+
+ /* Add diagnostic markers from markers structure */
+ for (i = 0; ; i++) {
+ double pp[3];
+ co cp;
+ if (markers[i].type == 0)
+ break;
+
+ if (markers[i].type == 1) { /* Src point - do luminance mapping */
+ dopartialmap1(s, pp, markers[i].pos);
+ } else {
+ pp[0] = markers[i].pos[0];
+ pp[1] = markers[i].pos[1];
+ pp[2] = markers[i].pos[2];
+ }
+ wrl->add_marker(wrl, pp, markers[i].col, 1.0);
+ }
+ }
+
+ /* --------------------------- */
+ /* Now computue our 3D mapping points from the near point mapping. */
+ for (i = 0; i < nnsm; i++) {
+ double cpexf; /* The effective compression or expansion factor */
+
+ 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");
+ }
+
+ /* Compute destination value which is a blend */
+ /* between the source value and the fully mapped destination value. */
+ icmBlend3(nsm[i].div, nsm[i].sv, nsm[i].dv, cpexf);
+
+#ifdef NEVER
+ printf("%s mapping:\n",nsm[i].vflag == 0 ? "Unclear" : nsm[i].vflag == 1 ? "Compression" : "Expansion");
+ printf("Src point = %f %f %f radius %f\n",nsm[i].sv[0], nsm[i].sv[1], nsm[i].sv[2], nsm[i].sr);
+ printf("Dst point = %f %f %f radius %f\n",nsm[i].dv[0], nsm[i].dv[1], nsm[i].dv[2], nsm[i].dr);
+ printf("Blended dst point = %f %f %f\n",nsm[i].div[0], nsm[i].div[1], nsm[i].div[2]);
+#endif /* NEVER */
+ /* Set the main gamut hull mapping point */
+ for (j = 0; j < 3; j++) {
+ gpnts[ngamp].p[j] = nsm[i].sv[j];
+ gpnts[ngamp].v[j] = nsm[i].div[j];
+ }
+ gpnts[ngamp++].w = 1.01; /* Main gamut surface mapping point */
+ /* (Use 1.01 as a marker value) */
+
+#ifdef USE_GAMKNF
+ /* Add sub surface mapping point if available */
+ if (nsm[i].vflag != 0) { /* Sub surface point is available */
+
+ /* Compute destination value which is a blend */
+ /* between the source value and the fully mapped destination value. */
+ icmBlend3(nsm[i].div2, nsm[i].sv2, nsm[i].dv2, cpexf);
+
+#ifdef NEVER
+ printf("Src2 point = %f %f %f radius %f\n",nsm[i].sv2[0], nsm[i].sv2[1], nsm[i].sv2[2], nsm[i].sr);
+ printf("Dst2 point = %f %f %f radius %f\n",nsm[i].dv2[0], nsm[i].dv2[1], nsm[i].dv2[2], nsm[i].dr);
+ printf("Blended dst2 point = %f %f %f\n",nsm[i].div2[0], nsm[i].div2[1], nsm[i].div2[2]);
+ printf("\n");
+#endif /* NEVER */
+ /* Set the sub-surface gamut hull mapping point */
+ for (j = 0; j < 3; j++) {
+ gpnts[ngamp].p[j] = nsm[i].sv2[j];
+ gpnts[ngamp].v[j] = nsm[i].div2[j];
+ }
+ gpnts[ngamp++].w = nsm[i].w2; /* Sub-suface mapping points */
+ }
+#endif /* USE_GAMKNF */
+ }
+
+ /* 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;
+ }
+ 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++) {
+ for (j = 0; j < 3; j++) {
+ if (gpnts[i].p[j] < (il[j]-1e-5) || gpnts[i].p[j] > (ih[j]+1e-5))
+ warning("gammap internal: input bounds has changed! %f <> %f <> %f",il[j],gpnts[i].p[j],ih[j]);
+ if (gpnts[i].v[j] < ol[j])
+ ol[j] = gpnts[i].v[j];
+ if (gpnts[i].v[j] > oh[j])
+ oh[j] = gpnts[i].v[j];
+ }
+ }
+
+ /* --------------------------- */
+
+#ifdef NEVER /* Dump out all the mapping points */
+ {
+ for (i = 0; i < ngamp; i++) {
+ printf("%d: %f %f %f -> %f %f %f\n",i,
+ gpnts[i].p[0], gpnts[i].p[1], gpnts[i].p[2],
+ gpnts[i].v[0], gpnts[i].v[1], gpnts[i].v[2]);
+ }
+ }
+#endif
+
+ /* Create the final gamut mapping rspl. */
+ /* [ 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)
+ printf("Creating rspl..\n");
+ for (j = 0; j < 3; j++) { /* Set resolution for all axes */
+ gres[j] = mapres;
+ avgdev[j] = GAMMAP_RSPLAVGDEV;
+ }
+ s->map = new_rspl(RSPL_NOFLAGS, 3, 3); /* Allocate 3D -> 3D */
+ if (s->map->fit_rspl_w(s->map, GAMMAP_RSPLFLAGS, gpnts, ngamp, il, ih, gres, ol, oh, smooth, avgdev, NULL)) {
+ if (verb)
+ fprintf(stderr,"Warning: Gamut mapping is non-monotonic - may not be very smooth !\n");
+ }
+ /* return the min and max of the input values valid in the grid */
+ s->map->get_in_range(s->map, s->imin, s->imax);
+
+#ifdef CHECK_NEARMAP
+ /* Check how accurate gamut shell mapping is against nsm */
+ /* (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; /* DE stats */
+
+ for (i = 0; i < nnsm; i++) {
+ double av[3];
+
+ /* Compute the mapping error */
+ dopartialmap2(s, av, nsm[i].sv); /* Just the rspl */
+
+ de = icmLabDE(nsm[i].div, av);
+ avgde += de;
+ if (de > maxde)
+ maxde = de;
+ }
+ printf("Gamut hull fit to guides: = avg %f, max %f\n",avgde/nnsm,maxde);
+ }
+#endif /* CHECK_NEARMAP */
+
+ /* If requested, enhance the saturation of the output values. */
+ if (gmi->satenh > 0.0) {
+ adjustsat cx; /* Adjustment context */
+
+ /* Compute what our source white and black points actually maps to */
+ s->domap(s, cx.wp, s_mt_wp);
+ s->domap(s, cx.bp, s_mt_bp);
+
+ cx.dst = d_gam;
+ cx.satenh = gmi->satenh;
+
+ /* Saturation enhance the output values */
+ s->map->re_set_rspl(
+ s->map, /* this */
+ 0, /* Combination of flags */
+ (void *)&cx, /* Opaque function context */
+ adjust_sat_func /* Function to set from */
+ );
+ }
+
+ /* Test the gamut white and black point mapping, and "fine tune" */
+ /* the mapping, to ensure an accurate transform of the white */
+ /* and black points to the destination colorspace. */
+ /* This compensates for any inacuracy introduced in the */
+ /* various rspl mappings. */
+ {
+ adjustwb cx; /* Adjustment context */
+ double a_wp[3]; /* actual white point */
+ double a_bp[3]; /* actual black point */
+
+ if (verb)
+ printf("Fine tuning white and black point mapping\n");
+
+ /* Check what the source white and black points actually maps to */
+ s->domap(s, a_wp, s_mt_wp);
+ s->domap(s, a_bp, s_mt_bp);
+
+#ifdef VERBOSE
+ if (verb) {
+ printf("White is %f %f %f, should be %f %f %f\n",
+ a_wp[0], a_wp[1], a_wp[2], d_mt_wp[0], d_mt_wp[1], d_mt_wp[2]);
+ printf("Black is %f %f %f, should be %f %f %f\n",
+ a_bp[0], a_bp[1], a_bp[2], d_mt_bp[0], d_mt_bp[1], d_mt_bp[2]);
+ }
+#endif /* VERBOSE */
+
+ /* Setup the fine tune transform */
+
+ /* We've decided not to fine tune the black point if we're */
+ /* bending to the destination black, as the bend is not */
+ /* followed perfectly (too sharp, or in conflict with */
+ /* the surface mapping ?) and we don't want to shift */
+ /* mid neutrals due to this. */
+ /* We do fine tune it if dst_kbp is set though, since */
+ /* we would like perfect K only out. */
+
+ /* Compute rotation/scale relative white point matrix */
+ icmVecRotMat(cx.mat, a_wp, a_bp, d_mt_wp, d_mt_bp); /* wp & bp */
+
+ /* Fine tune the 3D->3D mapping */
+ s->map->re_set_rspl(
+ s->map, /* this */
+ 0, /* Combination of flags */
+ (void *)&cx, /* Opaque function context */
+ adjust_wb_func /* Function to set from */
+ );
+
+#ifdef VERBOSE
+ if (verb) {
+ /* Check what the source white and black points actually maps to */
+ s->domap(s, a_wp, s_mt_wp);
+ s->domap(s, a_bp, s_mt_bp);
+
+ printf("After fine tuning:\n");
+ printf("White is %f %f %f, should be %f %f %f\n",
+ a_wp[0], a_wp[1], a_wp[2], d_mt_wp[0], d_mt_wp[1], d_mt_wp[2]);
+ printf("Black is %f %f %f, should be %f %f %f\n",
+ a_bp[0], a_bp[1], a_bp[2], d_mt_bp[0], d_mt_bp[1], d_mt_bp[2]);
+ }
+#endif /* VERBOSE */
+ }
+
+ if (wrl != NULL) {
+ int arerings = 0;
+ double cc[3] = { 0.7, 0.7, 0.7 };
+ double nc[3] = { 1.0, 0.4, 0.7 }; /* Pink for neighbors */
+ int nix = -1; /* Index of point to show neighbour */
+
+#ifdef SHOW_NEIGBORS
+#ifdef NEVER
+ /* Show all neighbours */
+ wrl->start_line_set(wrl, 0);
+ for (i = 0; i < nnsm; i++) {
+ for (j = 0; j < XNNB; j++) {
+ nearsmth *np = nsm[i].n[j]; /* Pointer to neighbor */
+
+ if (np == NULL)
+ break;
+
+ wrl->add_col_vertex(wrl, 0, nsm[i].sv, nc); /* Source value */
+ wrl->add_col_vertex(wrl, 0, np->sv, nc); /* Neighbpor value */
+ }
+ }
+ wrl->make_lines(wrl, 0, 2);
+#else
+ /* Show neighbours of points near source markers */
+ for (i = 0; ; i++) { /* Add diagnostic markers */
+ double pp[3];
+ co cp;
+ int ix, bix;
+ double bdist = 1e6;
+
+ if (markers[i].type == 0)
+ break;
+
+ if (markers[i].type != 1)
+ continue;
+
+ /* Rotate and map marker point the same as the src gamuts */
+ icmMul3By3x4(pp, s->grot, markers[i].pos);
+ cp.p[0] = pp[0]; /* L value */
+ s->grey->interp(s->grey, &cp);
+ pp[0] = cp.v[0];
+//printf("~1 looking for closest point to marker %d at %f %f %f\n",i,pp[0],pp[1],pp[2]);
+
+ /* Locate the nearest source point */
+ for (ix = 0; ix < nnsm; ix++) {
+ double dist = icmNorm33(pp, nsm[ix].sv);
+ if (dist < bdist) {
+ bdist = dist;
+ bix = ix;
+ }
+ }
+//printf("~1 closest src point ix %d at %f %f %f\n",bix,nsm[bix].sv[0],nsm[bix].sv[1],nsm[bix].sv[2]);
+//printf("~1 there are %d neighbours\n",nsm[bix].nnb);
+
+ wrl->start_line_set(wrl, 0);
+ for (j = 0; j < nsm[bix].nnb; j++) {
+ nearsmth *np = nsm[bix].n[j].n; /* Pointer to neighbor */
+
+ wrl->add_col_vertex(wrl, 0, nsm[bix].sv, nc); /* Source value */
+ wrl->add_col_vertex(wrl, 0, np->sv, nc); /* Neighbpor value */
+ }
+ wrl->make_lines(wrl, 0, 2);
+ }
+#endif
+#endif /* SHOW_NEIGBORS */
+
+ /* Add the source and dest gamut surfaces */
+#ifdef PLOT_SRC_GMT
+ wrl->make_gamut_surface_2(wrl, sil_gam, 0.6, 0, cc);
+#endif /* PLOT_SRC_GMT */
+#ifdef PLOT_DST_GMT
+ cc[0] = -1.0;
+ wrl->make_gamut_surface(wrl, d_gam, 0.2, cc);
+#endif /* PLOT_DST_GMT */
+#ifdef PLOT_DIGAM
+ if (nsm[0].dgam == NULL)
+ error("Need to #define PLOT_DIGAM in nearsmth.c!");
+ cc[0] = -1.0;
+ wrl->make_gamut_surface(wrl, nsm[0].dgam, 0.2, cc);
+#endif /* PLOT_DIGAM */
+#ifdef PLOT_SRC_CUSPS
+ wrl->add_cusps(wrl, sil_gam, 0.6, NULL);
+#endif /* PLOT_SRC_CUSPS */
+#ifdef PLOT_DST_CUSPS
+ wrl->add_cusps(wrl, d_gam, 0.3, NULL);
+#endif /* PLOT_DST_CUSPS */
+
+#ifdef PLOT_TRANSSRC_CUSPS
+ /* Add transformed source cusp markers */
+ {
+ int i;
+ double cusps[6][3];
+ 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 */
+ };
+
+ if (sc_gam->getcusps(sc_gam, cusps) == 0) {
+
+ for (i = 0; i < 6; i++) {
+ double val[3];
+
+ s->domap(s, val, cusps[i]);
+ wrl->add_marker(wrl, val, ccolors[i], 2.5);
+ }
+ }
+ }
+#endif
+
+#if defined(SHOW_MAP_VECTORS) || defined(SHOW_SUB_SURF) || defined(SHOW_ACTUAL_VECTORS) || defined(SHOW_ACTUAL_VEC_DIFF)
+ /* Start of guide vector plot */
+ wrl->start_line_set(wrl, 0);
+
+ for (i = 0; i < nnsm; i++) {
+ double cpexf; /* The effective compression or expansion factor */
+ 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 lgrey[3] = { 0.8, 0.8, 0.8 };
+ double purp[3] = { 0.6, 0.0, 1.0 };
+ double blue[3] = { 0.2, 0.2, 1.0 };
+ double *ccc;
+ double mdst[3];
+
+#if defined(SHOW_ACTUAL_VECTORS) || defined(SHOW_ACTUAL_VEC_DIFF)
+# ifdef SHOW_ACTUAL_VECTORS
+ wrl->add_col_vertex(wrl, 0, nsm[i].sv, yellow);
+# else /* SHOW_ACTUAL_VEC_DIFF */
+ wrl->add_col_vertex(wrl, 0, nsm[i].div, yellow);
+# endif
+ dopartialmap2(s, mdst, nsm[i].sv);
+ wrl->add_col_vertex(wrl, 0, mdst, red);
+
+#else
+# ifdef SHOW_MAP_VECTORS
+ ccc = yellow;
+
+ if (nsm[i].gflag == 0)
+ ccc = green; /* Mark "no clear direction" vectors in green->red */
+# ifdef SHOW_CUSPMAP
+ wrl->add_col_vertex(wrl, 0, nsm[i].csv, ccc); /* Cusp mapped source value */
+# else
+ wrl->add_col_vertex(wrl, 0, nsm[i].sv, ccc); /* Source value */
+# endif
+ wrl->add_col_vertex(wrl, 0, nsm[i].div, red); /* Blended destination value */
+# endif /* SHOW_MAP_VECTORS */
+
+# ifdef SHOW_SUB_SURF
+ if (nsm[i].vflag != 0) { /* Sub surface point is available */
+
+ wrl->add_col_vertex(wrl, 0, nsm[i].sv2, lgrey); /* Subs-surf Source value */
+ wrl->add_col_vertex(wrl, 0, nsm[i].div2, purp); /* Blended destination value */
+ }
+# endif /* SHOW_SUB_SURF */
+#endif /* !SHOW_ACTUAL_VECTORS */
+ }
+ wrl->make_lines(wrl, 0, 2); /* Guide vectors */
+#endif /* Show vectors */
+
+#ifdef SHOW_VECTOR_INDEXES
+ for (i = 0; i < nnsm; i++) {
+ double cream[3] = { 0.7, 0.7, 0.5 };
+ char buf[100];
+ sprintf(buf, "%d", i);
+ wrl->add_text(wrl, buf, nsm[i].sv, cream, 0.5);
+ }
+#endif /* SHOW_VECTOR_INDEXES */
+
+ /* add the locus from locus.ts file */
+ if (locus != NULL) {
+ int table, npoints;
+ char *fnames[3] = { "LAB_L", "LAB_A", "LAB_B" };
+ int ix[3];
+ double v0[3], v1[3];
+ double rgb[3];
+
+ /* Each table holds a separate locus */
+ for (table = 0; table < locus->ntables; table++) {
+
+ if ((npoints = locus->t[table].nsets) <= 0)
+ error("No sets of data in diagnostic locus");
+
+ for (j = 0; j < 3; j++) {
+ if ((ix[j] = locus->find_field(locus, 0, fnames[j])) < 0)
+ error ("Locus file doesn't contain field %s",fnames[j]);
+ if (locus->t[table].ftype[ix[j]] != r_t)
+ error ("Field %s is wrong type",fnames[j]);
+ }
+
+ /* Source locus */
+ rgb[0] = 1.0;
+ rgb[1] = 0.5;
+ rgb[2] = 0.5;
+ for (i = 0; i < npoints; i++) {
+ co cp;
+
+ for (j = 0; j < 3; j++)
+ v1[j] = *((double *)locus->t[table].fdata[i][ix[j]]);
+
+ /* Rotate and locus verticies the same as the src gamuts */
+ dopartialmap1(s, v1, v1);
+ if (i > 0 )
+ wrl->add_cone(wrl, v0, v1, rgb, 0.5);
+ icmAry2Ary(v0,v1);
+ }
+
+ /* Gamut mapped locus */
+ rgb[0] = 1.0;
+ rgb[1] = 1.0;
+ rgb[2] = 1.0;
+ for (i = 0; i < npoints; i++) {
+ co cp;
+
+ for (j = 0; j < 3; j++)
+ v1[j] = *((double *)locus->t[table].fdata[i][ix[j]]);
+
+ s->domap(s, v1, v1);
+ if (i > 0 )
+ wrl->add_cone(wrl, v0, v1, rgb, 0.5);
+ icmAry2Ary(v0,v1);
+ }
+ }
+
+ locus->del(locus);
+ locus = NULL;
+ }
+
+ /* Add any ring mapping diagnostics */
+ for (i = 0; ; i++) {
+ if (rings[i].type == 0)
+ break;
+
+ if (rings[i].type == 2)
+ continue;
+
+ if (rings[i].type == 1) {
+ double pconst;
+ double cpoint[3];
+ double mat[3][4]; /* translate to our plane */
+ double imat[3][4]; /* translate from our plane */
+ double s1[3], s0[3], t1[3];
+ int j;
+ double maxa, mina;
+ double maxb, minb;
+
+ if (arerings == 0) {
+ arerings = 1;
+ wrl->start_line_set(wrl, 1); /* Source ring */
+ wrl->start_line_set(wrl, 2); /* Destination ring */
+ }
+
+ if (icmNormalize3(rings[i].pnorm, rings[i].pnorm, 1.0))
+ error("Ring %d diagnostic plane normal failed",i);
+
+ pconst = -icmDot3(rings[i].ppoint, rings[i].pnorm);
+
+ /* Locate intersection of source neautral axis and plane */
+ if (icmVecPlaneIsect(cpoint, pconst, rings[i].pnorm, s_cs_wp, s_cs_bp))
+ error("Ring %d diagnostic center point intersection failed",i);
+
+ /* Compute the rotation and translation between */
+ /* a plane in ab and the plane we are using */
+ s0[0] = s0[1] = s0[2] = 0.0;
+ s1[0] = 1.0, s1[1] = s1[2] = 0.0;
+ t1[0] = cpoint[0] + rings[i].pnorm[0];
+ t1[1] = cpoint[1] + rings[i].pnorm[1];
+ t1[2] = cpoint[2] + rings[i].pnorm[2];
+ icmVecRotMat(mat, s1, s0, t1, cpoint);
+ icmVecRotMat(imat, t1, cpoint, s1, s0);
+
+ /* Do a min/max of a circle of vectors so as to */
+ /* establish an offset to the centroid for this slice */
+ maxa = maxb = -1e60;
+ mina = minb = 1e60;
+ for (j = 0; j < 20; j++) {
+ double ang = 2 * 3.1415926 * j/(20 - 1.0);
+ double vec[3], isect[3];
+ double pp[3];
+ co cp;
+ int k;
+
+ vec[0] = 0.0;
+ vec[1] = sin(ang);
+ vec[2] = cos(ang);
+ icmMul3By3x4(vec, mat, vec);
+
+ /* Intersect it with the source gamut */
+ if (si_gam->vector_isect(si_gam, vec, cpoint, isect,
+ NULL, NULL, NULL, NULL, NULL) == 0) {
+ continue;
+ }
+
+ /* Translate back to plane */
+ icmMul3By3x4(pp, imat, isect);
+
+ if (pp[1] > maxa)
+ maxa = pp[1];
+ if (pp[1] < mina)
+ mina = pp[1];
+ if (pp[2] > maxb)
+ maxb = pp[2];
+ if (pp[2] < minb)
+ minb = pp[2];
+ }
+ /* Move center to centroid of min/max box */
+ t1[0] = 0.0;
+ t1[1] = (maxa + mina) * 0.5;
+ t1[2] = (maxb + minb) * 0.5;
+ if (t1[1] < -200.0 || t1[1] > 200.0
+ || t1[2] < -200.0 || t1[2] > 200.0)
+ error("Failed to locate centroid of slice");
+ icmMul3By3x4(cpoint, mat, t1);
+
+//printf("~1 ring centroid point = %f %f %f\n", cpoint[0],cpoint[1],cpoint[2]);
+
+ /* Recompute the rotation and translation between */
+ /* a plane in ab and the plane we are using */
+ s0[0] = s0[1] = s0[2] = 0.0;
+ s1[0] = 1.0, s1[1] = s1[2] = 0.0;
+ t1[0] = cpoint[0] + rings[i].pnorm[0];
+ t1[1] = cpoint[1] + rings[i].pnorm[1];
+ t1[2] = cpoint[2] + rings[i].pnorm[2];
+ icmVecRotMat(mat, s1, s0, t1, cpoint);
+ icmVecRotMat(imat, t1, cpoint, s1, s0);
+
+//printf("~1 generating %d ring verts\n",rings[i].nverts);
+ /* Create a circle of vectors in the plane from the center */
+ /* point, to intersect with the source gamut surface. */
+ /* (Duplicate start and end vertex) */
+ for (j = 0; j <= rings[i].nverts; j++) {
+ double ang = 2 * 3.1415926 * j/((double) rings[i].nverts);
+ double vec[3], isect[3];
+ double pp[3];
+ co cp;
+ int k;
+
+ vec[0] = 0.0;
+ vec[1] = sin(ang);
+ vec[2] = cos(ang);
+ icmMul3By3x4(vec, mat, vec);
+
+ /* Intersect it with the source gamut */
+ 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;
+ }
+
+//printf("~1 vec %d = %f %f %f\n",j,isect[0],isect[1],isect[2]);
+
+ /* Scale them to the ratio */
+ for (k = 0; k < 3; k++)
+ vec[k] = isect[k] * rings[i].rad + (1.0 - rings[i].rad) * cpoint[k];
+
+//printf("~1 rad vec %d = %f %f %f\n",j,vec[0],vec[1],vec[2]);
+
+ /* Transform them into rotated and scaled destination space */
+ dopartialmap1(s, vec, vec);
+//printf("~1 trans vec %d = %f %f %f\n",j,vec[0],vec[1],vec[2]);
+
+ /* Add to plot */
+ wrl->add_col_vertex(wrl, 1, vec, rings[i].scol);
+//printf("~1 src vec %d = %f %f %f\n",j,vec[0],vec[1],vec[2]);
+
+ /* Gamut map and add to plot */
+ s->domap(s, vec, vec);
+//printf("~1 dst vec %d = %f %f %f\n",j,vec[0],vec[1],vec[2]);
+ wrl->add_col_vertex(wrl, 2, vec, rings[i].dcol);
+ }
+ wrl->make_last_vertex(wrl, 1); /* Source ring */
+ wrl->make_last_vertex(wrl, 2); /* Destination ring */
+ }
+ if (arerings) {
+ wrl->make_lines(wrl, 1, 1000000); /* Source ring */
+ wrl->make_lines(wrl, 2, 1000000); /* Destination ring */
+ }
+ }
+
+ wrl->del(wrl); /* Write and delete */
+ wrl = NULL;
+ }
+
+#ifdef PLOT_3DKNEES
+ /* Plot one graph per 3D gamut boundary mapping point */
+ for (j = 0; j < p3dk_ix; j++) {
+ double xx[XRES];
+ double yy[XRES];
+
+ printf("Vector %f %f %f -> %f %f %f\n", p3dk_locus[j].v0[0], p3dk_locus[j].v0[1], p3dk_locus[j].v0[2], p3dk_locus[j].v1[0], p3dk_locus[j].v1[1], p3dk_locus[j].v1[2]);
+
+ for (i = 0; i < XRES; i++) {
+ double v;
+ co cp; /* Conversion point */
+ v = (i/(double)(XRES-1.0));
+ cp.p[0] = p3dk_locus[j].v0[0] + v * (p3dk_locus[j].v1[0] - p3dk_locus[j].v0[0]);
+ cp.p[1] = p3dk_locus[j].v0[1] + v * (p3dk_locus[j].v1[1] - p3dk_locus[j].v0[1]);
+ cp.p[2] = p3dk_locus[j].v0[2] + v * (p3dk_locus[j].v1[2] - p3dk_locus[j].v0[2]);
+ xx[i] = sqrt(cp.p[1] * cp.p[1] + cp.p[2] * cp.p[2]);
+ s->map->interp(s->map, &cp);
+ yy[i] = sqrt(cp.v[1] * cp.v[1] + cp.v[2] * cp.v[2]);
+ }
+ do_plot(xx,yy,NULL,NULL,XRES);
+ }
+ free(p3dk_locus);
+#endif /* PLOT_3DKNEES */
+
+
+ free(gpnts);
+ free_nearsmth(nsm, nnsm);
+
+ } else if (diagname != NULL && verb) {
+ printf("Warning: Won't create '%s' because there is no 3D gamut mapping\n",diagname);
+ }
+
+#ifdef PLOT_GAMUTS
+ scl_gam->write_vrml(scl_gam, "src.wrl", 1, 0);
+ sil_gam->write_vrml(sil_gam, "img.wrl", 1, 0);
+ d_gam->write_vrml(d_gam, "dst.wrl", 1, 0);
+ sc_gam->write_trans_vrml(sc_gam, "gmsrc.wrl", 1, 0, map_trans, s);
+#endif
+
+ if (sil_gam != scl_gam)
+ sil_gam->del(sil_gam);
+ scl_gam->del(scl_gam);
+
+ return s;
+}
+
+#ifdef PLOT_GAMUTS
+
+/* Debug */
+static void map_trans(void *cntx, double out[3], double in[3]) {
+ gammap *map = (gammap *)cntx;
+
+ map->domap(map, out, in);
+}
+
+#endif
+
+/* Object methods */
+static void del_gammap(
+gammap *s
+) {
+ if (s->grey != NULL)
+ s->grey->del(s->grey);
+ if (s->igrey != NULL)
+ s->igrey->del(s->igrey);
+ if (s->map != NULL)
+ s->map->del(s->map);
+
+ free(s);
+}
+
+/* Apply the gamut mapping to the given color value */
+static void domap(
+gammap *s,
+double *out,
+double *in
+) {
+ double rin[3];
+ co cp;
+
+ if (s->dbg) printf("domap: got input %f %f %f\n",in[0],in[1],in[2]);
+ icmMul3By3x4(rin, s->grot, in); /* Rotate */
+
+ if (s->dbg) printf("domap: after rotate %f %f %f\n",rin[0],rin[1],rin[2]);
+ cp.p[0] = rin[0];
+ s->grey->interp(s->grey, &cp); /* L map */
+
+ if (s->dbg) printf("domap: after L map %f %f %f\n",cp.v[0],rin[1],rin[2]);
+
+ /* If there is a 3D->3D mapping */
+ if (s->map != NULL) {
+ int e;
+
+ /* Clip out of range a, b proportionately */
+ if (rin[1] < s->imin[1] || rin[1] > s->imax[1]
+ || rin[2] < s->imin[2] || rin[2] > s->imax[2]) {
+ double as = 1.0, bs = 1.0;
+ if (rin[1] < s->imin[1])
+ as = s->imin[1]/rin[1];
+ else if (rin[1] > s->imax[1])
+ as = s->imax[1]/rin[1];
+ if (rin[2] < s->imin[2])
+ bs = s->imin[2]/rin[2];
+ else if (rin[2] > s->imax[2])
+ bs = s->imax[2]/rin[2];
+ if (bs < as)
+ as = bs;
+ rin[1] *= as;
+ rin[2] *= as;
+ }
+
+ cp.p[0] = cp.v[0]; /* 3D map */
+ cp.p[1] = rin[1];
+ cp.p[2] = rin[2];
+ s->map->interp(s->map, &cp);
+
+ for (e = 0; e < s->map->fdi; e++)
+ out[e] = cp.v[e];
+
+ if (s->dbg) printf("domap: after 3D map %s\n\n",icmPdv(s->map->fdi, out));
+ } else {
+ out[0] = cp.v[0];
+ out[1] = rin[1];
+ out[2] = rin[2];
+ }
+}
+
+/* Apply the matrix and grey mapping to the given color value */
+static void dopartialmap1(
+gammap *s,
+double *out,
+double *in
+) {
+ double rin[3];
+ co cp;
+
+ icmMul3By3x4(rin, s->grot, in); /* Rotate */
+ cp.p[0] = rin[0];
+ s->grey->interp(s->grey, &cp); /* L map */
+ out[0] = cp.v[0];
+ out[1] = rin[1];
+ out[2] = rin[2];
+}
+
+/* Apply just the rspl mapping to the given color value */
+/* (ie. to a color already rotated and L mapped) */
+static void dopartialmap2(
+gammap *s,
+double *out,
+double *in
+) {
+ co cp;
+
+ /* If there is a 3D->3D mapping */
+ if (s->map != NULL) {
+ int e;
+
+ icmCpy3(cp.p, in);
+
+ /* Clip out of range a, b proportionately */
+ if (cp.p[1] < s->imin[1] || cp.p[1] > s->imax[1]
+ || cp.p[2] < s->imin[2] || cp.p[2] > s->imax[2]) {
+ double as = 1.0, bs = 1.0;
+ if (cp.p[1] < s->imin[1])
+ as = s->imin[1]/cp.p[1];
+ else if (cp.p[1] > s->imax[1])
+ as = s->imax[1]/cp.p[1];
+ if (cp.p[2] < s->imin[2])
+ bs = s->imin[2]/cp.p[2];
+ else if (cp.p[2] > s->imax[2])
+ bs = s->imax[2]/cp.p[2];
+ if (bs < as)
+ as = bs;
+ cp.p[1] *= as;
+ cp.p[2] *= as;
+ }
+
+ s->map->interp(s->map, &cp);
+
+ icmCpy3(out, cp.v);
+ } else {
+ icmCpy3(out, in);
+ }
+}
+
+/* Function to pass to rspl to invert grey curve */
+static void inv_grey_func(
+ void *cntx,
+ double *out,
+ double *in
+) {
+ rspl *fwd = (rspl *)cntx;
+ int nsoln; /* Number of solutions found */
+ co pp[2]; /* Room for all the solutions found */
+
+ pp[0].p[0] =
+ pp[0].v[0] = in[0];
+
+ nsoln = fwd->rev_interp(
+ fwd,
+ RSPL_NEARCLIP, /* Clip to nearest (faster than vector) */
+ 2, /* Maximum number of solutions allowed for */
+ NULL, /* No auxiliary input targets */
+ NULL, /* Clip vector direction and length */
+ pp); /* Input and output values */
+
+ nsoln &= RSPL_NOSOLNS; /* Get number of solutions */
+
+ if (nsoln != 1)
+ error("gammap: Unexpected failure to find reverse solution for grey axis lookup");
+
+ out[0] = pp[0].p[0];
+}
+
+/* Function to pass to rspl to alter output values, */
+/* to enhance the saturation. */
+static void
+adjust_sat_func(
+ void *pp, /* adjustsat structure */
+ double *out, /* output value to be adjusted */
+ double *in /* corresponding input value */
+) {
+ adjustsat *p = (adjustsat *)pp;
+ double cp[3]; /* Center point */
+ double rr, t1[3], p1;
+ double t2[3], p2;
+
+ /* Locate center point on the white/black axis corresponding to this color */
+ cp[0] = out[0];
+ rr = (out[0] - p->bp[0])/(p->wp[0] - p->bp[0]); /* Relative location on the white/black axis */
+ cp[1] = p->bp[1] + rr * (p->wp[1] - p->bp[1]);
+ cp[2] = p->bp[2] + rr * (p->wp[2] - p->bp[2]);
+
+ /* Locate the point on the destination gamut surface in the direction */
+ /* from the center point to the point being processed. */
+ if (p->dst->vector_isect(p->dst, cp, out, t2, t1, &p2, &p1, NULL, NULL) != 0) {
+
+ if (p1 > 1.0) { /* If this point is within gamut */
+ double ep1, bf;
+
+//printf("\n");
+//printf("~1 cp %f %f %f input %f %f %f\n",cp[0],cp[1],cp[2], out[0], out[1], out[2]);
+//printf("~1 min %f %f %f mint %f\n",t2[0],t2[1],t2[2],p2);
+//printf("~1 max %f %f %f maxt %f\n",t1[0],t1[1],t1[2],p1);
+
+ p1 = 1.0/p1; /* Position of out from cp to t1 */
+
+#ifdef NEVER
+ /* Enhanced parameter value */
+ ep1 = (p1 + p->satenh * p1)/(1.0 + p->satenh * p1);
+ /* Make blend between linear p1 and enhanced p1, */
+ /* to reduce effects on near neutrals. */
+ p1 = (1.0 - p1) * p1 + p1 * ep1;
+#else
+ /* Compute Enhanced p1 */
+ ep1 = (p1 + p->satenh * p1)/(1.0 + p->satenh * p1);
+
+ /* Make blend factor between linear p1 and enhanced p1, */
+ /* to reduce effects on near neutrals. */
+ {
+ double pp = 4.0; /* Sets where the 50% transition is */
+ double g = 2.0; /* Sets rate of transition */
+ double sec, vv = p1;
+
+ vv = vv/(pp - pp * vv + 1.0);
+
+ vv *= 2.0;
+ sec = floor(vv);
+ if (((int)sec) & 1)
+ g = -g; /* Alternate action in each section */
+ vv -= sec;
+ if (g >= 0.0) {
+ vv = vv/(g - g * vv + 1.0);
+ } else {
+ vv = (vv - g * vv)/(1.0 - g * vv);
+ }
+ vv += sec;
+ vv *= 0.5;
+
+ bf = (vv + pp * vv)/(1.0 + pp * vv);
+ }
+ /* Do the blend */
+ p1 = (1.0 - bf) * p1 + bf * ep1;
+#endif
+ /* Compute enhanced values position */
+ out[0] = cp[0] + (t1[0] - cp[0]) * p1;
+ out[1] = cp[1] + (t1[1] - cp[1]) * p1;
+ out[2] = cp[2] + (t1[2] - cp[2]) * p1;
+//printf("~1 output %f %f %f, param %f\n",out[0],out[1],out[2],p1);
+ }
+ }
+}
+
+/* Function to pass to rspl to re-set output values, */
+/* to adjust the white and black points */
+static void
+adjust_wb_func(
+ void *pp, /* adjustwb structure */
+ double *out, /* output value to be adjusted */
+ double *in /* corresponding input value */
+) {
+ adjustwb *p = (adjustwb *)pp;
+
+ /* Do a linear mapping from swp -> dwp and sbp -> dbp, */
+ /* to compute the adjusted value. */
+ icmMul3By3x4(out, p->mat, out);
+}
+
+
+/* Create a new gamut that the the given gamut transformed by the */
+/* gamut mappings rotation and grey curve mapping. Return NULL on error. */
+static gamut *parttransgamut(gammap *s, gamut *src) {
+ gamut *dst;
+ double cusps[6][3];
+ double wp[3], bp[3], kp[3];
+ double p[3];
+ int i;
+
+ if ((dst = new_gamut(src->getsres(src), src->getisjab(src), src->getisrast(src))) == NULL)
+ return NULL;
+
+ dst->setnofilt(dst);
+
+ /* Translate all the surface nodes */
+ for (i = 0;;) {
+ if ((i = src->getrawvert(src, p, i)) < 0)
+ break;
+
+ dopartialmap1(s, p, p);
+ dst->expand(dst, p);
+ }
+ /* Translate cusps */
+ if (src->getcusps(src, cusps) == 0) {
+ dst->setcusps(dst, 0, NULL);
+ for (i = 0; i < 6; i++) {
+ dopartialmap1(s, p, cusps[i]);
+ dst->setcusps(dst, 1, p);
+ }
+ dst->setcusps(dst, 2, NULL);
+ }
+ /* Translate white and black points */
+ if (src->getwb(src, wp, bp, kp, NULL, NULL, NULL) == 0) {
+ dopartialmap1(s, wp, wp);
+ dopartialmap1(s, bp, bp);
+ dopartialmap1(s, kp, kp);
+ dst->setwb(dst, wp, bp, kp);
+ }
+
+ return dst;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-01 00:00:06 +0000