summaryrefslogtreecommitdiff
path: root/xicc/cam02.c
diff options
context:
space:
mode:
Diffstat (limited to 'xicc/cam02.c')
-rw-r--r--xicc/cam02.c1279
1 files changed, 1279 insertions, 0 deletions
diff --git a/xicc/cam02.c b/xicc/cam02.c
new file mode 100644
index 0000000..1ab6abc
--- /dev/null
+++ b/xicc/cam02.c
@@ -0,0 +1,1279 @@
+
+/*
+ * cam02
+ *
+ * Color Appearance Model, based on
+ * CIECAM02, "The CIECAM02 Color Appearance Model"
+ * by Nathan Moroney, Mark D. Fairchild, Robert W.G. Hunt, Changjun Li,
+ * M. Ronnier Luo and Todd Newman, IS&T/SID Tenth Color Imaging
+ * Conference, with the addition of the Viewing Flare
+ * model described on page 487 of "Digital Color Management",
+ * by Edward Giorgianni and Thomas Madden, and the
+ * Helmholtz-Kohlraush effect, using the equation from
+ * the Bradford-Hunt 96C model as detailed in Mark Fairchild's
+ * book "Color Appearance Models".
+ * The Slight modification to the Hunt-Pointer-Estevez matrix
+ * recommended by Kuo, Zeis and Lai in IS&T/SID 14th Color Imaging
+ * Conference "Robust CIECAM02 Implementation and Numerical
+ * Experiment within an ICC Workflow", page 215, together with
+ * their matrix formulation of inversion has been adopted.
+ *
+ * In addition the algorithm has unique extensions to allow
+ * it to operate reasonably over an unbounded domain.
+ *
+ * Author: Graeme W. Gill
+ * Date: 17/1/2004
+ * Version: 1.00
+ *
+ * This file is based on cam97s3.c by Graeme Gill.
+ *
+ * Copyright 2004 - 2011 Graeme W. Gill
+ * Please refer to COPYRIGHT file for details.
+ * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
+ * see the License.txt file for licencing details.
+ */
+
+
+/* Note that XYZ values are normalised to 1.0 consistent */
+/* with the ICC convention (not 100.0 as assumed by the CIECAM spec.) */
+/* Note that all whites are assumed to be normalised (ie. Y = 1.0) */
+
+/*
+ Various additions and changes have been made to allow the CAM
+ conversions to and from an unbounded range of XYZ and Jab values,
+ in a (somewhat) geometrically consistent maner. This is because
+ such values arise in the process of gamut mapping, and in scanning through
+ the grid of PCS values needed to fill in the A2B table of an
+ ICC profile. Such values have no correlation to a real color
+ value, but none the less need to be handled without causing an
+ exception, in a geometrically consistent and reversible
+ fashion.
+
+ The changes are:
+
+ The Hunt-Pointer-Estevez matrix has been increased in precission by
+ 1 digit [Kuo, Zeise & Lai, IS&T/SID 14th Color Imaging Conference pp. 215-219]
+
+ The sharpened cone response matrix third row has been changed from
+ [0.0030, 0.0136, 0.9834] to [0.0000, 0.0000, 1.0000]
+ [Susstrunk & Brill, IS&T/SID 14th Color Imaging Conference Late Breaking News sublement]
+
+ To prevent wild Jab values as the rgb' values approach zero, a region close to each
+ primary ground plane is compressed. This expands the region of reasonable
+ behaviour to encompass XYZ values that are found in practice (ie. ProPhotoRGB).
+
+ To prevent excessive curvature of hue in high saturation blue regions
+ due to the non-linearity exagerating the difference between
+ r & b values, a heuristic is introduced that blends the r & b
+ values, reducing this effect. This degrades the agreement
+ with the reference CIECAM implementation by about 1DE in this region,
+ but greatly improves the behaviour in mapping and clipping from
+ highly saturated blues (ie. ProPhotoRGB).
+
+ The post-adaptation non-linearity response has had a straight
+ line negative linear extension added.
+ The positive region has a tangential linear extension added, so
+ that the range of values is not limited.
+ An adaptive threshold for the low level linear extension,
+ to avoid coordinates blowing up when one of the cone responses
+ is large, while the others become negative.
+
+ Re-arrange ss equation into separated effects,
+ k1, k2, k3, and then limit these to avoid divide by zero
+ and sign change errors in fwd and bwd converson, as well
+ as avoiding the blue saturation doubling back for low
+ luminance levels.
+
+ To avoid chroma and hue angle information being lost when the
+ J value used to scale the chroma is 0, and to ensure
+ that J = 0, a,b != 0 values have a valid mapping into
+ XYZ space, the J value used to multiply Chroma, is limited
+ to be equivalent to not less than A == 0.1.
+
+ The Helmholtz-Kohlraush effect is crafted to have resonable
+ effects over the range of J from 0 to 100, and have more
+ moderate effects outside this range.
+
+*/
+
+#include <copyright.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "icc.h"
+#include "xcam.h"
+#include "cam02.h"
+#include "numlib.h"
+
+#define ENABLE_COMPR /* [Def] Enable XYZ compression */
+#define ENABLE_BLUE_ANGLE_FIX /* [Def] Limit maximum blue angle */
+#define ENABLE_DDL /* [Def] Enable k1,k2,k3 overall ss limit values (seems to be the best scheme) */
+#undef ENABLE_SS /* [Undef] Disable overall ss limit values (not the scheme used) */
+
+#undef ENTRACE /* [Undef] Enable internal value runtime tracing if s->trace != 0 */
+#undef DOTRACE /* [Undef] Trace anyway (ie. set s->trace = 1) */
+#undef DIAG1 /* [Undef] Print internal value diagnostics for conditions setup */
+#undef DIAG2 /* [Undef] Print internal value diagnostics for each conversion */
+#undef TRACKMINMAX /* [Undef] Track min/max DD & SS limits (run with locus cam02test) */
+#undef DISABLE_MATRIX /* Debug - disable matrix & non-lin, wire XYZ to rgba */
+#undef DISABLE_SCR /* Debug - disable Sharpened Cone Response matrix */
+#undef DISABLE_HPE /* Debug - disable just Hunt-Pointer_Estevez matrix */
+#undef DISABLE_NONLIN /* Debug - wire rgbp to rgba */
+#undef DISABLE_TTD /* Debug - disable ttd vector 'tilt' */
+#undef DISABLE_HHKR /* Debug - disable Helmholtz-Kohlraush */
+
+#ifdef ENABLE_COMPR
+# define BC_WHMINY 0.3 /* [0.3] Compression direction minimum Y value */
+# define BC_RANGE_R 0.01 /* [0.01] Set compression range as prop. of distance to neutral - red */
+# define BC_RANGE_G 0.05 /* [0.05] Set compression range as prop. of distance to neutral - green*/
+# define BC_RANGE_B 0.10 /* [0.10] Set compression range as prop. of distance to neutral - blue */
+# define BC_MAXRANGE 0.13 /* [0.13] Maximum compression range */
+# define BC_LIMIT 0.7 /* [0.7] Correction limit (abs. rgbp distance shift) */
+#endif /* ENABLE_COMPR */
+
+#ifdef ENABLE_BLUE_ANGLE_FIX
+# define BLUE_BL_MAX 0.9 /* [0.9] Sets ultimate blue angle, higher = less angle */
+# define BLUE_BL_POW 3.5 /* [3.5] Sets rate at which blue angle is limited */
+#endif
+
+#define NLDLIMIT 0.00001 /* [0.00001] Non-linearity minimum lower crossover to straight line */
+#define NLDICEPT -0.18 /* [-0.18] Input intercept of straight line with 0.1 output */
+
+#define NLULIMIT 1e5 /* Non-linearity upper crossover to straight line */
+
+#ifdef ENABLE_SS /* [Undef] */
+# define SSLLIMIT 0.22 /* Overall ab scale lower limit */
+# define SSULIMIT 580.0 /* Overall ab scale upper limit */
+#endif /* ENABLE_SS */
+
+#define SYMETRICJ /* [Undef] Use symetric power about zero, else straigt line -ve */
+
+#define DDLLIMIT 0.55 /* [0.55] ab component -k3:k2 ratio limit (must be < 1.0) */
+//#define DDULIMIT 0.9993 /* [0.9993] ab component k3:k1 ratio limit (must be < 1.0) */
+#define DDULIMIT 0.34 /* [0.34] ab component k3:k1 ratio limit (must be < 1.0) */
+#define SSMINcJ 0.005 /* [0.005] ab scale cJ minimum value */
+#define JLIMIT 0.005 /* [0.005] J encoding cutover point straight line (0 - 1.0 range) */
+#define HKLIMIT 0.7 /* [0.7] Maximum Helmholtz-Kohlraush lift */
+
+#ifdef TRACKMINMAX
+double minss = 1e60;
+double maxss = -1e60;
+double minlrat = 0.0;
+double maxurat = 0.0;
+#define noslots 103
+double slotsd[noslots];
+double slotsu[noslots];
+double minj = 1e38, maxj = -1e38;
+#endif /* TRACKMINMAX */
+
+#if defined(ENTRACE) || defined(DOTRACE)
+#define TRACE(xxxx) if (s->trace) printf xxxx ;
+#else
+#define TRACE(xxxx)
+#endif
+
+static void cam_free(cam02 *s);
+static int set_view(struct _cam02 *s, ViewingCondition Ev, double Wxyz[3],
+ double La, double Yb, double Lv, double Yf, double Fxyz[3],
+ int hk);
+static int XYZ_to_cam(struct _cam02 *s, double *Jab, double *xyz);
+static int cam_to_XYZ(struct _cam02 *s, double *xyz, double *Jab);
+
+static double spow(double val, double pp) {
+ if (val < 0.0)
+ return -pow(-val, pp);
+ else
+ return pow(val, pp);
+}
+
+/* Create a cam02 conversion object, with default viewing conditions */
+cam02 *new_cam02(void) {
+ cam02 *s;
+// double D50[3] = { 0.9642, 1.0000, 0.8249 };
+
+ if ((s = (cam02 *)calloc(1, sizeof(cam02))) == NULL) {
+ fprintf(stderr,"cam02: malloc failed allocating object\n");
+ exit(-1);
+ }
+
+ /* Initialise methods */
+ s->del = cam_free;
+ s->set_view = set_view;
+ s->XYZ_to_cam = XYZ_to_cam;
+ s->cam_to_XYZ = cam_to_XYZ;
+
+ /* Set default range handling limits */
+ s->nldlimit = NLDLIMIT;
+ s->nldicept = NLDICEPT;
+ s->nlulimit = NLULIMIT;
+ s->ddllimit = DDLLIMIT;
+ s->ddulimit = DDULIMIT;
+ s->ssmincj = SSMINcJ;
+ s->jlimit = JLIMIT;
+ s->hklimit = 1.0 / HKLIMIT;
+
+ /* Set a default viewing condition ?? */
+ /* set_view(s, vc_average, D50, 33, 0.2, 0.0, 0.0, D50, 0); */
+
+#ifdef DOTRACE
+ s->trace = 1;
+#endif
+
+#ifdef TRACKMINMAX
+ {
+ int i;
+ for (i = 0; i < noslots; i++) {
+ slotsd[i] = 0.0;
+ slotsu[i] = 0.0;
+ }
+ }
+#endif /* TRACKMINMAX */
+
+ return s;
+}
+
+static void cam_free(cam02 *s) {
+
+#ifdef TRACKMINMAX
+ {
+ int i;
+ for (i = 0; i < noslots; i++) {
+ printf("Slot %d = %f, %f\n",i,slotsd[i], slotsu[i]);
+ }
+ printf("minj = %f, maxj = %f\n",minj,maxj);
+
+ printf("minss = %f\n",minss);
+ printf("maxss = %f\n",maxss);
+ printf("minlrat = %f\n",minlrat);
+ printf("maxurat = %f\n",maxurat);
+ }
+#endif /* TRACKMINMAX */
+
+ if (s != NULL)
+ free(s);
+}
+
+/* Return value is always 0 */
+static int set_view(
+cam02 *s,
+ViewingCondition Ev, /* Enumerated Viewing Condition */
+double Wxyz[3], /* Reference/Adapted White XYZ (Y range 0.0 .. 1.0) */
+double La, /* Adapting/Surround Luminance cd/m^2 */
+double Yb, /* Relative Luminance of Background to reference white (range 0.0 .. 1.0) */
+double Lv, /* Luminance of white in the Viewing/Scene/Image field (cd/m^2) */
+ /* Ignored if Ev is set to other than vc_none */
+double Yf, /* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
+double Fxyz[3], /* The Flare white coordinates (typically the Ambient color) */
+int hk /* Flag, NZ to use Helmholtz-Kohlraush effect */
+) {
+ double tt, t1, t2;
+ double tm[3][3];
+ int i;
+
+ if (Ev == vc_none) {
+ /* Compute the internal parameters by interpolation */
+ int i;
+ double r, bf;
+ /* Dark, dim, average, above average */
+ double t_C[4] = { 0.525, 0.59, 0.69, 1.0 };
+ double t_Nc[4] = { 0.800, 0.95, 1.00, 1.0 };
+ double t_F[4] = { 0.800, 0.90, 1.00, 1.0 };
+
+ if (La < 1e-10) /* Hmm. */
+ La = 1e-10;
+ r = La/Lv;
+ if (r < 0.0)
+ r = 0.0;
+ else if (r > 1.0)
+ r = 1.0;
+
+ if (r < 0.1) { /* Dark to Dim */
+ i = 0;
+ bf = r/0.1;
+ } else if (r < 0.2) { /* Dim to Average */
+ i = 1;
+ bf = (r - 0.1)/0.1;
+ } else { /* Average to above average */
+ i = 2;
+ bf = (r - 0.2)/0.8;
+ }
+ s->C = t_C[i] * (1.0 - bf) + t_C[i+1] * bf;
+ s->Nc = t_Nc[i] * (1.0 - bf) + t_Nc[i+1] * bf;
+ s->F = t_F[i] * (1.0 - bf) + t_F[i+1] * bf;
+ } else {
+ /* Compute the internal parameters by category */
+ switch(Ev) {
+ case vc_dark:
+ s->C = 0.525;
+ s->Nc = 0.8;
+ s->F = 0.8;
+ break;
+ case vc_dim:
+ s->C = 0.59;
+ s->Nc = 0.95;
+ s->F = 0.9;
+ break;
+ case vc_cut_sheet:
+ s->C = 0.41;
+ s->Nc = 0.8;
+ s->F = 0.8;
+ break;
+ default: /* average */
+ s->C = 0.69;
+ s->Nc = 1.0;
+ s->F = 1.0;
+ break;
+ }
+ }
+
+ /* Transfer parameters to the object */
+ s->Ev = Ev;
+ s->Wxyz[0] = Wxyz[0];
+ s->Wxyz[1] = Wxyz[1];
+ s->Wxyz[2] = Wxyz[2];
+ s->Yb = Yb > 0.005 ? Yb : 0.005; /* Set minimum to avoid divide by 0.0 */
+ s->La = La;
+ s->Yf = Yf;
+ s->Fxyz[0] = Fxyz[0];
+ s->Fxyz[1] = Fxyz[1];
+ s->Fxyz[2] = Fxyz[2];
+ s->hk = hk;
+
+ /* The rgba vectors */
+ s->Va[0] = 1.0;
+ s->Va[1] = -12.0/11.0;
+ s->Va[2] = 1.0/11.0;
+
+ s->Vb[0] = 1.0/9.0;
+ s->Vb[1] = 1.0/9.0;
+ s->Vb[2] = -2.0/9.0;
+
+ s->VttA[0] = 2.0;
+ s->VttA[1] = 1.0;
+ s->VttA[2] = 1.0/20.0;
+
+ /* (Not used) */
+ s->Vttd[0] = 1.0;
+ s->Vttd[1] = 1.0;
+ s->Vttd[2] = 21.0/20.0;
+
+ /* Vttd in terms of the VttA, Va and Vb vectors */
+ s->dcomp[0] = 1.0;
+ s->dcomp[1] = -11.0/23.0;
+ s->dcomp[2] = -108.0/23.0;
+
+ /* Compute values that only change with viewing parameters */
+
+ /* Figure out the Flare contribution to the flareless XYZ input */
+ tt = s->Yf * s->Wxyz[1]/s->Fxyz[1];
+ s->Fsxyz[0] = tt * s->Fxyz[0];
+ s->Fsxyz[1] = tt * s->Fxyz[1];
+ s->Fsxyz[2] = tt * s->Fxyz[2];
+
+ /* Rescale so that the sum of the flare and the input doesn't exceed white */
+ s->Fsc = s->Wxyz[1]/(s->Fsxyz[1] + s->Wxyz[1]);
+ s->Fsxyz[0] *= s->Fsc;
+ s->Fsxyz[1] *= s->Fsc;
+ s->Fsxyz[2] *= s->Fsc;
+ s->Fisc = 1.0/s->Fsc;
+
+#ifndef DISABLE_SCR
+ /* Sharpened cone response white values */
+ s->rgbW[0] = 0.7328 * s->Wxyz[0] + 0.4296 * s->Wxyz[1] - 0.1624 * s->Wxyz[2];
+ s->rgbW[1] = -0.7036 * s->Wxyz[0] + 1.6975 * s->Wxyz[1] + 0.0061 * s->Wxyz[2];
+ s->rgbW[2] = 0.0000 * s->Wxyz[0] + 0.0000 * s->Wxyz[1] + 1.0000 * s->Wxyz[2];
+#else
+ s->rgbW[0] = s->Wxyz[0];
+ s->rgbW[1] = s->Wxyz[1];
+ s->rgbW[2] = s->Wxyz[2];
+#endif
+
+ /* Degree of chromatic adaptation */
+ s->D = s->F * (1.0 - exp((-s->La - 42.0)/92.0)/3.6);
+
+ /* Precompute Chromatic transform values */
+ s->Drgb[0] = s->D * (s->Wxyz[1]/s->rgbW[0]) + 1.0 - s->D;
+ s->Drgb[1] = s->D * (s->Wxyz[1]/s->rgbW[1]) + 1.0 - s->D;
+ s->Drgb[2] = s->D * (s->Wxyz[1]/s->rgbW[2]) + 1.0 - s->D;
+
+ /* Chromaticaly transformed white value */
+ s->rgbcW[0] = s->Drgb[0] * s->rgbW[0];
+ s->rgbcW[1] = s->Drgb[1] * s->rgbW[1];
+ s->rgbcW[2] = s->Drgb[2] * s->rgbW[2];
+
+#ifndef DISABLE_HPE
+ /* Transform from spectrally sharpened, to Hunt-Pointer_Estevez cone space */
+ s->rgbpW[0] = 0.7409744840453773 * s->rgbcW[0]
+ + 0.2180245944753982 * s->rgbcW[1]
+ + 0.0410009214792244 * s->rgbcW[2];
+ s->rgbpW[1] = 0.2853532916858801 * s->rgbcW[0]
+ + 0.6242015741188157 * s->rgbcW[1]
+ + 0.0904451341953042 * s->rgbcW[2];
+ s->rgbpW[2] = -0.0096276087384294 * s->rgbcW[0]
+ - 0.0056980312161134 * s->rgbcW[1]
+ + 1.0153256399545427 * s->rgbcW[2];
+#else
+ s->rgbpW[0] = s->rgbcW[0];
+ s->rgbpW[1] = s->rgbcW[1];
+ s->rgbpW[2] = s->rgbcW[2];
+#endif /* DISABLE_HPE */
+
+#ifndef DISABLE_SCR
+ /* Create combined cone and chromatic transform matrix: */
+ /* Spectrally sharpened cone responses matrix */
+ s->cc[0][0] = 0.7328; s->cc[0][1] = 0.4296; s->cc[0][2] = -0.1624;
+ s->cc[1][0] = -0.7036; s->cc[1][1] = 1.6975; s->cc[1][2] = 0.0061;
+ s->cc[2][0] = 0.0000; s->cc[2][1] = 0.0000; s->cc[2][2] = 1.0000;
+#else
+ s->cc[0][0] = 1.0; s->cc[0][1] = 0.0; s->cc[0][2] = 0.0;
+ s->cc[1][0] = 0.0; s->cc[1][1] = 1.0; s->cc[1][2] = 0.0;
+ s->cc[2][0] = 0.0; s->cc[2][1] = 0.0; s->cc[2][2] = 1.0;
+#endif
+
+ /* Chromaticaly transformed sample value */
+ icmSetUnity3x3(tm);
+ tm[0][0] = s->Drgb[0];
+ tm[1][1] = s->Drgb[1];
+ tm[2][2] = s->Drgb[2];
+ icmMul3x3(s->cc, tm);
+
+#ifndef DISABLE_HPE
+ /* Transform from spectrally sharpened, to Hunt-Pointer_Estevez cone space */
+ tm[0][0] = 0.7409744840453773;
+ tm[0][1] = 0.2180245944753982;
+ tm[0][2] = 0.0410009214792244;
+ tm[1][0] = 0.2853532916858801;
+ tm[1][1] = 0.6242015741188157;
+ tm[1][2] = 0.0904451341953042;
+ tm[2][0] = -0.0096276087384294;
+ tm[2][1] = -0.0056980312161134;
+ tm[2][2] = 1.0153256399545427;
+ icmMul3x3(s->cc, tm);
+#endif /* DISABLE_HPE */
+
+ /* Create inverse combined cone and chromatic transform matrix: */
+ icmInverse3x3(s->icc, s->cc); /* Hmm. we assume it cannot fail */
+
+#ifdef ENABLE_COMPR
+ /* Compression ranges */
+ s->crange[0] = BC_RANGE_R;
+ s->crange[1] = BC_RANGE_G;
+ s->crange[2] = BC_RANGE_B;
+#endif /* ENABLE_COMPR */
+
+ /* Background induction factor */
+ s->n = s->Yb/ s->Wxyz[1];
+ s->nn = pow(1.64 - pow(0.29, s->n), 0.73); /* Pre computed value */
+
+ /* Lightness contrast factor ?? */
+ {
+ double k;
+
+ k = 1.0 / (5.0 * s->La + 1.0);
+ s->Fl = 0.2 * pow(k , 4.0) * 5.0 * s->La
+ + 0.1 * pow(1.0 - pow(k , 4.0) , 2.0) * pow(5.0 * s->La , 1.0/3.0);
+ }
+
+ /* Background and Chromatic brightness induction factors */
+ s->Nbb = 0.725 * pow(1.0/s->n, 0.2);
+ s->Ncb = s->Nbb;
+
+ /* Base exponential nonlinearity */
+ s->z = 1.48 + pow(s->n , 0.5);
+
+ /* Post-adapted cone response of white */
+ tt = pow(s->Fl * s->rgbpW[0], 0.42);
+ s->rgbaW[0] = 400.0 * tt / (tt + 27.13) + 0.1;
+ tt = pow(s->Fl * s->rgbpW[1], 0.42);
+ s->rgbaW[1] = 400.0 * tt / (tt + 27.13) + 0.1;
+ tt = pow(s->Fl * s->rgbpW[2], 0.42);
+ s->rgbaW[2] = 400.0 * tt / (tt + 27.13) + 0.1;
+
+ /* Achromatic response of white */
+ s->Aw = (s->VttA[0] * s->rgbaW[0]
+ + s->VttA[1] * s->rgbaW[1]
+ + s->VttA[2] * s->rgbaW[2] - 0.305) * s->Nbb;
+
+ /* Non-linearity lower crossover output value */
+ tt = pow(s->Fl * s->nldlimit, 0.42);
+ s->nldxval = 400.0 * tt / (tt + 27.13) + 0.1;
+
+ /* Non-linearity lower crossover slope from lower crossover */
+ /* to intercept with 0.1 output */
+ s->nldxslope = (s->nldxval - 0.1)/(s->nldlimit - s->nldicept);
+
+ /* Non-linearity upper crossover value */
+ tt = pow(s->Fl * s->nlulimit, 0.42);
+ s->nluxval = 400.0 * tt / (tt + 27.13) + 0.1;
+
+ /* Non-linearity upper crossover slope, set to asymtope */
+ t1 = s->Fl * s->nlulimit;
+ t2 = pow(t1, 0.42) + 27.13;
+ s->nluxslope = 0.42 * s->Fl * 400.0 * 27.13 / (pow(t1,0.58) * t2 * t2);
+
+
+ /* Limited A value at J = JLIMIT */
+ s->lA = pow(s->jlimit, 1.0/(s->C * s->z)) * s->Aw;
+
+#ifdef DIAG1
+ printf("Scene parameters:\n");
+ printf("Viewing condition Ev = %d\n",s->Ev);
+ printf("Ref white Wxyz = %f %f %f\n", s->Wxyz[0], s->Wxyz[1], s->Wxyz[2]);
+ printf("Relative liminance of background Yb = %f\n", s->Yb);
+ printf("Adapting liminance La = %f\n", s->La);
+ printf("Flare Yf = %f\n", s->Yf);
+ printf("Flare color Fxyz = %f %f %f\n", s->Fxyz[0], s->Fxyz[1], s->Fxyz[2]);
+
+ printf("Internal parameters:\n");
+ printf("Surround Impact C = %f\n", s->C);
+ printf("Chromatic Induction Nc = %f\n", s->Nc);
+ printf("Adaptation Degree F = %f\n", s->F);
+
+ printf("Pre-computed values\n");
+ printf("Sharpened cone white rgbW = %f %f %f\n", s->rgbW[0], s->rgbW[1], s->rgbW[2]);
+ printf("Degree of chromatic adaptation D = %f\n", s->D);
+ printf("Chromatic transform values Drgb = %f %f %f\n", s->Drgb[0], s->Drgb[1], s->Drgb[2]);
+ printf("Chromatically transformed white rgbcW = %f %f %f\n", s->rgbcW[0], s->rgbcW[1], s->rgbcW[2]);
+ printf("Hunter-P-E cone response white rgbpW = %f %f %f\n", s->rgbpW[0], s->rgbpW[1], s->rgbpW[2]);
+ printf("Background induction factor n = %f\n", s->n);
+ printf(" nn = %f\n", s->nn);
+ printf("Lightness contrast factor Fl = %f\n", s->Fl);
+ printf("Background brightness induction factor Nbb = %f\n", s->Nbb);
+ printf("Chromatic brightness induction factor Ncb = %f\n", s->Ncb);
+ printf("Base exponential nonlinearity z = %f\n", s->z);
+ printf("Post adapted cone response white rgbaW = %f %f %f\n", s->rgbaW[0], s->rgbaW[1], s->rgbaW[2]);
+ printf("Achromatic response of white Aw = %f\n", s->Aw);
+ printf("\n");
+#endif
+ return 0;
+}
+
+/* Conversions. Return values are always 0 */
+static int XYZ_to_cam(
+struct _cam02 *s,
+double Jab[3],
+double XYZ[3]
+) {
+ int i;
+ double xyz[3], rgbp[3], rgba[3];
+ double a, b, ja, jb, J, JJ, C, h, e, A, ss;
+ double ttA, rS, cJ, tt;
+ double k1, k2, k3;
+ int clip = 0;
+
+#ifdef DIAG2 /* Incase in == out */
+ double XYZi[3];
+
+ XYZi[0] = XYZ[0];
+ XYZi[1] = XYZ[1];
+ XYZi[2] = XYZ[2];
+#endif
+
+ TRACE(("\nForward conversion:\n"))
+ TRACE(("XYZ = %f %f %f\n",XYZ[0], XYZ[1], XYZ[2]))
+
+#ifdef DISABLE_MATRIX
+
+ rgba[0] = XYZ[0];
+ rgba[1] = XYZ[1];
+ rgba[2] = XYZ[2];
+
+#else /* !DISABLE_MATRIX */
+
+ /* Add in flare */
+ xyz[0] = s->Fsc * XYZ[0] + s->Fsxyz[0];
+ xyz[1] = s->Fsc * XYZ[1] + s->Fsxyz[1];
+ xyz[2] = s->Fsc * XYZ[2] + s->Fsxyz[2];
+
+ TRACE(("XYZ inc flare = %f %f %f\n",xyz[0], xyz[1], xyz[2]))
+
+ /* Spectrally sharpened cone responses, */
+ /* Chromaticaly transformed sample value, */
+ /* Transform from spectrally sharpened, to Hunt-Pointer_Estevez cone space. */
+ icmMulBy3x3(rgbp, s->cc, xyz);
+
+ TRACE(("rgbp = %f %f %f\n", rgbp[0], rgbp[1], rgbp[2]))
+
+#ifdef ENABLE_COMPR
+ /* Try and prevent crazy out of cam02 gamut behaviour, by compressing */
+ /* the rgbp so as to prevent it becoming less than zero. */
+ {
+ double tt; /* Temporary */
+ double wrgb[3]; /* White target */
+
+ /* Make white target white point with same Y value */
+ tt = xyz[1] > BC_WHMINY ? xyz[1] : BC_WHMINY; /* Limit to minimum Y */
+ icmScale3(wrgb, s->rgbpW, tt/s->Wxyz[1]); /* White target at same Y */
+ TRACE(("wrgb %f %f %f\n", wrgb[0], wrgb[1], wrgb[2]))
+
+ /* Compress r,g,b in turn */
+ for (i = 0; i < 3; i++) {
+ double cv; /* Compression value */
+ double ctv; /* Compression target value */
+ double cd; /* Compression displacement needed */
+ double cvec[3]; /* Normalized correction vector */
+ double isec[3]; /* Intersection with plane */
+ double offs; /* Offset of point from orgin*/
+ double range; /* Threshold to start compression */
+ double asym; /* Compression asymtope */
+
+ /* Compute compression direction as vector towards white */
+ /* (We did try correcting in a blend of limit plane normal and white, */
+ /* but compressing towards white seems to be the best.) */
+ icmSub3(cvec, wrgb, rgbp); /* Direction of white target */
+
+ TRACE(("rgbp %f %f %f\n", rgbp[0], rgbp[1], rgbp[2]))
+ TRACE(("cvec %f %f %f\n", cvec[0], cvec[1], cvec[2]))
+
+ if (cvec[i] < 1e-9) { /* compression direction can't correct this coord */
+ TRACE(("Intersection with limit plane failed\n"))
+ continue;
+ }
+
+ /* Scale compression vector to make it move a unit in normal direction */
+ icmScale3(cvec, cvec, 1.0/cvec[i]); /* Normalized vector to white */
+ TRACE(("cvec %f %f %f\n", cvec[0], cvec[1], cvec[2]))
+
+ /* Compute intersection of correction direction with this limit plane */
+ /* (This corresponds with finding displacement of rgbp by cvec */
+ /* such that the current coord value = 0) */
+ icmScale3(isec, cvec, -rgbp[i]); /* (since cvec[i] == 1.0) */
+ icmAdd3(isec, isec, rgbp);
+ TRACE(("isec %f %f %f\n", isec[0], isec[1], isec[2]))
+
+ /* Compute distance from intersection to origin */
+ offs = pow(icmNorm3(isec), 0.85);
+
+ range = s->crange[i] * offs; /* Scale range by distance to origin */
+ if (range > BC_MAXRANGE) /* so that it tapers down as we approach it */
+ range = BC_MAXRANGE; /* and limit maximum */
+
+ /* Aiming above plane when far from origin, */
+ /* but below plane at the origin, so that black isn't affected. */
+ asym = range - 0.2 * (range + (0.01 * s->crange[i]));
+
+ ctv = cv = rgbp[i]; /* Distance above/below limit plane */
+
+ TRACE(("ch %d, offs %f, range %f asym %f, cv %f\n",i, offs,range,asym,cv))
+ if (cv < (range - 1e-12)) { /* Need to compress */
+ double aa, bb;
+ aa = 1.0/(range - cv);
+ bb = 1.0/(range - asym);
+ ctv = range - 1.0/(aa + bb);
+
+ cd = ctv - cv; /* Displacement needed */
+ if (cd > BC_LIMIT)
+ cd = BC_LIMIT;
+ TRACE(("ch %d cd = %f, scaled cd %f\n",i,cd,cd))
+
+ /* Apply correction */
+ icmScale3(cvec, cvec, cd); /* Compression vector */
+ icmAdd3(rgbp, rgbp, cvec); /* Compress by displacement */
+ TRACE(("rgbp after comp. = %f %f %f\n",rgbp[0], rgbp[1], rgbp[2]))
+ }
+ }
+ }
+#endif /* ENABLE_COMPR */
+
+#ifdef ENABLE_BLUE_ANGLE_FIX
+ ss = rgbp[0] + rgbp[1] + rgbp[2];
+ if (ss < 1e-9) {
+ ss = 0.0;
+ } else {
+ ss = (rgbp[2]/ss - 1.0/3.0) * 3.0/2.0;
+ if (ss > 0.0)
+ ss = BLUE_BL_MAX * pow(ss, BLUE_BL_POW);
+ }
+ if (ss < 0.0)
+ ss = 0.0;
+ else if (ss > 1.0)
+ ss = 1.0;
+ tt = 0.5 * (rgbp[0] + rgbp[1]);
+ rgbp[0] = ss * tt + (1.0 - ss) * rgbp[0];
+ rgbp[1] = ss * tt + (1.0 - ss) * rgbp[1];
+ TRACE(("rgbp after blue fix ss = %f, rgbp = %f %f %f\n",ss,rgbp[0], rgbp[1], rgbp[2]))
+#endif
+
+#ifdef DISABLE_NONLIN
+ for (i = 0; i < 3; i++) {
+ rgba[i] = 400.0/27.13 * rgbp[i];
+ }
+#else /* !DISABLE_NONLIN */
+ /* Post-adapted cone response of sample. */
+ /* rgba[] has a minimum value of 0.1 for XYZ[] = 0 and no flare. */
+ /* We add a negative linear region, plus a linear segment at */
+ /* the end of the +ve conversion to allow numerical handling of a */
+ /* very wide range of values. */
+
+ for (i = 0; i < 3; i++) {
+ if (rgbp[i] < s->nldlimit) {
+ rgba[i] = s->nldxval + s->nldxslope * (rgbp[i] - s->nldlimit);
+ } else {
+ if (rgbp[i] <= s->nlulimit) {
+ tt = pow(s->Fl * rgbp[i], 0.42);
+ rgba[i] = 400.0 * tt / (tt + 27.13) + 0.1;
+ } else {
+ rgba[i] = s->nluxval + s->nluxslope * (rgbp[i] - s->nlulimit);
+ }
+ }
+ }
+#endif /* !DISABLE_NONLIN */
+
+//tt = 0.5 * (rgba[0] + rgba[1]);
+//rgba[0] = (rgba[0] - ss * tt)/(1.0 - ss);
+//rgba[1] = (rgba[1] - ss * tt)/(1.0 - ss);
+
+#endif /* !DISABLE_MATRIX */
+
+ /* Note that the minimum values of rgba[] for XYZ = 0 is 0.1, */
+ /* hence magic 0.305 below comes from the following weighting of rgba[], */
+ /* to base A at 0.0 */
+ /* Deal with the core rgb to A, S & h conversion: */
+
+ TRACE(("rgba = %f %f %f\n", rgba[0], rgba[1], rgba[2]))
+
+ /* Preliminary Acromatic response */
+ ttA = s->VttA[0] * rgba[0] + s->VttA[1] * rgba[1] + s->VttA[2] * rgba[2];
+
+ /* Achromatic response */
+ A = (ttA - 0.305) * s->Nbb;
+
+ /* Preliminary red-green & yellow-blue opponent dimensions */
+ /* a, b & ttd form an (almost) orthogonal coordinate set. */
+ /* ttA is in a different direction */
+ a = s->Va[0] * rgba[0] + s->Va[1] * rgba[1] + s->Va[2] * rgba[2];
+ b = s->Vb[0] * rgba[0] + s->Vb[1] * rgba[1] + s->Vb[2] * rgba[2];
+
+ /* restricted Saturation to stop divide by zero */
+ /* (The exact value isn't important because the numerator dominates as a,b aproach 0 */
+ rS = sqrt(a * a + b * b);
+ if (rS < DBL_EPSILON)
+ rS = DBL_EPSILON;
+ TRACE(("ttA = %f, a = %f, b = %f, rS = %f, A = %f\n", ttA,a,b,rS,A))
+
+ /* Lightness J, Derived directly from Acromatic response. */
+ /* Cuttover to a straight line segment when J < 0.005, */
+#ifndef SYMETRICJ /* Cut to a straight line */
+ if (A >= s->lA) {
+ J = pow(A/s->Aw, s->C * s->z); /* J/100 - keep Sign */
+ } else {
+ J = s->jlimit/s->lA * A; /* Straight line */
+ TRACE(("limited Acromatic to straight line\n"))
+ }
+#else /* Symetric */
+ if (A >= 0.0) {
+ J = pow(A/s->Aw, s->C * s->z); /* J/100 - keep Sign */
+ } else {
+ J = -pow(-A/s->Aw, s->C * s->z); /* J/100 - keep Sign */
+ TRACE(("symetric Acromatic\n"))
+ }
+#endif
+
+ /* Constraied (+ve, non-zero) J */
+ if (A > 0.0) {
+ cJ = pow(A/s->Aw, s->C * s->z);
+ if (cJ < s->ssmincj)
+ cJ = s->ssmincj;
+ } else {
+ cJ = s->ssmincj;
+ }
+
+ TRACE(("J = %f, cJ = %f\n",J,cJ))
+
+ /* Final hue angle */
+ h = (180.0/DBL_PI) * atan2(b,a);
+ h = (h < 0.0) ? h + 360.0 : h;
+
+ /* Eccentricity factor */
+ e = (12500.0/13.0 * s->Nc * s->Ncb * (cos(h * DBL_PI/180.0 + 2.0) + 3.8));
+
+ /* ab scale components */
+ k1 = pow(s->nn, 1.0/0.9) * e * pow(cJ, 1.0/1.8)/pow(rS, 1.0/9.0);
+ k2 = pow(cJ, 1.0/(s->C * s->z)) * s->Aw/s->Nbb + 0.305;
+ k3 = s->dcomp[1] * a + s->dcomp[2] * b;
+
+ TRACE(("Raw k1 = %f, k2 = %f, k3 = %f, raw ss = %f\n",k1, k2, k3, pow(k1/(k2 + k3), 0.9)))
+
+#ifdef TRACKMINMAX
+ {
+ int sno;
+ double lrat, urat;
+
+ ss = pow(k1/(k2 + k3), 0.9);
+
+ if (ss < minss)
+ minss = ss;
+ if (ss > maxss)
+ maxss = ss;
+
+ lrat = -k3/k2;
+ urat = k3 * pow(ss, 10.0/9.0) / k1;
+ if (lrat > minlrat)
+ minlrat = lrat;
+ if (urat > maxurat)
+ maxurat = urat;
+
+ /* Record distribution of ss min/max vs. J for */
+ /* regions outside a,b == 0 */
+
+ sno = (int)(J * 100.0 + 0.5);
+ if (sno < 0) {
+ sno = 101;
+ if (J < minj)
+ minj = J;
+ } else if (sno > 100) {
+ sno = 102;
+ if (J > maxj)
+ maxj = J;
+ }
+ if (slotsd[sno] < lrat)
+ slotsd[sno] = lrat;
+ if (slotsu[sno] < urat)
+ slotsu[sno] = urat;
+ }
+#endif /* TRACKMINMAX */
+
+#ifdef ENABLE_DDL
+
+ /* Limit ratio of k3 to k2 to stop zero or -ve ss */
+ if (k3 < -k2 * s->ddllimit) {
+ k3 = -k2 * s->ddllimit;
+ TRACE(("k3 limited to %f due to k3:k2 ratio, ss = %f\n",k3,pow(k1/(k2 + k3), 0.9)))
+ clip = 1;
+ }
+
+ /* See if there is going to be a problem in bwd, and limit k3 if there is */
+ if (k3 > (k2 * s->ddulimit/(1.0 - s->ddulimit))) {
+ k3 = (k2 * s->ddulimit/(1.0 - s->ddulimit));
+ TRACE(("k3 limited to %f to allow for bk3:bk1 bwd limit\n",k3))
+ clip = 1;
+ }
+
+#endif /* !ENABLE_DDL */
+
+#ifdef DISABLE_TTD
+
+ ss = pow((k1/k2), 0.9);
+
+#else /* !TRACKMINMAX */
+
+ /* Compute the ab scale factor */
+ ss = pow(k1/(k2 + k3), 0.9);
+
+#endif /* !ENABLE_DDL */
+
+#ifdef ENABLE_SS
+ if (ss < SSLLIMIT)
+ ss = SSLLIMIT;
+ else if (ss > SSULIMIT)
+ ss = SSULIMIT;
+#endif /* ENABLE_SS */
+
+#ifdef NEVER
+// -------------------------------------------
+ /* Show ss components */
+ if (s->retss) {
+ Jab[0] = 1.0;
+ if (clip)
+ Jab[0] = 0.0;
+ Jab[1] = (log10(ss) - +1.5)/(2.0 - +1.5);
+ Jab[2] = (log10(ss) - +1.0)/(2.5 - +1.0);
+
+ for (i = 0; i < 3; i++) {
+ if (Jab[i] < 0.0)
+ Jab[i] = 0.0;
+ else if (Jab[i] > 1.0)
+ Jab[i] = 1.0;
+ }
+ return 0;
+ }
+// -------------------------------------------
+#endif /* NEVER */
+
+ ja = a * ss;
+ jb = b * ss;
+
+ /* Chroma - always +ve, used just for HHKR */
+ C = sqrt(ja * ja + jb * jb);
+
+ TRACE(("ss = %f, A = %f, J = %f, C = %f, h = %f\n",ss,A,J,C,h))
+
+ JJ = J;
+#ifndef DISABLE_HHKR
+ /* Helmholtz-Kohlraush effect */
+ if (s->hk && J < 1.0) {
+ double kk = C/300.0 * sin(DBL_PI * fabs(0.5 * (h - 90.0))/180.0);
+ if (kk > 1e-6) /* Limit kk to a reasonable range */
+ kk = 1.0/(s->hklimit + 1.0/kk);
+ JJ = J + (1.0 - (J > 0.0 ? J : 0.0)) * kk;
+ TRACE(("JJ = %f from J = %f, kk = %f\n",JJ,J,kk))
+ }
+#endif /* DISABLE_HHKR */
+
+ JJ *= 100.0; /* Scale J */
+
+ /* Compute Jab value */
+ Jab[0] = JJ;
+ Jab[1] = ja;
+ Jab[2] = jb;
+
+#ifdef NEVER /* Brightness/Colorfulness */
+ {
+ double M, Q;
+
+ Q = (4.0/s->C) * sqrt(JJ/100.0) * (s->Aw + 4.0) * pow(s->Fl, 0.25);
+ M = C * pow(s->Fl, 0.25);
+
+ printf("Lightness = %f, Chroma = %f\n",J,C);
+ printf("Brightness = %f, Colorfulness = %f\n",M,Q);
+ }
+#endif /* NEVER */
+
+ TRACE(("Jab %f %f %f\n",Jab[0], Jab[1], Jab[2]))
+ TRACE(("\n"))
+
+#ifdef DIAG2
+ printf("Processing XYZ->Jab:\n");
+ printf("XYZ = %f %f %f\n", XYZi[0], XYZi[1], XYZi[2]);
+ printf("Including flare XYZ = %f %f %f\n", xyz[0], xyz[1], xyz[2]);
+ printf("Huntpace rgbpP-E cone space rgbp = %f %f %f\n", rgbp[0], rgbp[1], rgbp[2]);
+ printf("Post adapted cone response rgba = %f %f %f\n", rgba[0], rgba[1], rgba[2]);
+ printf("Prelim red green a = %f, b = %f\n", a, b);
+ printf("Hue angle h = %f\n", h);
+ printf("Eccentricity factor e = %f\n", e);
+ printf("Achromatic response A = %f\n", A);
+ printf("Lightness J = %f, H.K. Lightness = %f\n", J * 100, JJ);
+ printf("Prelim Saturation ss = %f\n", ss);
+ printf("Chroma C = %f\n", C);
+ printf("Jab = %f %f %f\n", Jab[0], Jab[1], Jab[2]);
+ printf("\n");
+#endif
+ return 0;
+}
+
+static int cam_to_XYZ(
+struct _cam02 *s,
+double XYZ[3],
+double Jab[3]
+) {
+ int i;
+ double xyz[3], rgbp[3], rgba[3];
+ double a, b, ja, jb, J, JJ, C, rC, h, e, A, ss;
+ double tt, cJ, ttA;
+ double k1, k2, k3; /* (k1 & k3 are different to the fwd k1 & k3) */
+
+#ifdef DIAG2 /* Incase in == out */
+ double Jabi[3];
+
+ Jabi[0] = Jab[0];
+ Jabi[1] = Jab[1];
+ Jabi[2] = Jab[2];
+
+#endif
+
+ TRACE(("\nReverse conversion:\n"))
+ TRACE(("Jab %f %f %f\n",Jab[0], Jab[1], Jab[2]))
+
+ JJ = Jab[0] * 0.01; /* J/100 */
+ ja = Jab[1];
+ jb = Jab[2];
+
+ /* Convert Jab to core A, S & h values: */
+
+ /* Compute hue angle */
+ h = (180.0/DBL_PI) * atan2(jb, ja);
+ h = (h < 0.0) ? h + 360.0 : h;
+
+ /* Compute chroma value */
+ C = sqrt(ja * ja + jb * jb); /* Must be Always +ve, Can be NZ even if J == 0 */
+
+ /* Preliminary Restricted chroma, always +ve and NZ */
+ /* (The exact value isn't important because the numerator dominates as a,b aproach 0) */
+ rC = C;
+ if (rC < DBL_EPSILON)
+ rC = DBL_EPSILON;
+
+ J = JJ;
+#ifndef DISABLE_HHKR
+ /* Undo Helmholtz-Kohlraush effect */
+ if (s->hk && J < 1.0) {
+ double kk = C/300.0 * sin(DBL_PI * fabs(0.5 * (h - 90.0))/180.0);
+ if (kk > 1e-6) /* Limit kk to a reasonable range */
+ kk = 1.0/(s->hklimit + 1.0/kk);
+ J = (JJ - kk)/(1.0 - kk);
+ if (J < 0.0)
+ J = JJ - kk;
+ TRACE(("J = %f from JJ = %f, kk = %f\n",J,JJ,kk))
+ }
+#endif /* DISABLE_HHKR */
+
+ /* Achromatic response */
+#ifndef SYMETRICJ /* Cut to a straight line */
+ if (J >= s->jlimit) {
+ A = pow(J, 1.0/(s->C * s->z)) * s->Aw;
+ } else { /* In the straight line segment */
+ A = s->lA/s->jlimit * J;
+ TRACE(("Undo Acromatic straight line\n"))
+ }
+#else /* Symetric */
+ if (J >= 0.0) {
+ A = pow(J, 1.0/(s->C * s->z)) * s->Aw;
+ } else { /* In the straight line segment */
+ A = -pow(-J, 1.0/(s->C * s->z)) * s->Aw;
+ TRACE(("Undo symetric Acromatic\n"))
+ }
+#endif
+
+ /* Preliminary Acromatic response +ve */
+ ttA = (A/s->Nbb)+0.305;
+
+ if (A > 0.0) {
+ cJ = pow(A/s->Aw, s->C * s->z);
+ if (cJ < s->ssmincj)
+ cJ = s->ssmincj;
+ } else {
+ cJ = s->ssmincj;
+ }
+ TRACE(("C = %f, A = %f from J = %f, cJ = %f\n",C, A,J,cJ))
+
+ /* Eccentricity factor */
+ e = (12500.0/13.0 * s->Nc * s->Ncb * (cos(h * DBL_PI/180.0 + 2.0) + 3.8));
+
+ /* ab scale components */
+ k1 = pow(s->nn, 1.0/0.9) * e * pow(cJ, 1.0/1.8)/pow(rC, 1.0/9.0);
+ k2 = pow(cJ, 1.0/(s->C * s->z)) * s->Aw/s->Nbb + 0.305;
+ k3 = s->dcomp[1] * ja + s->dcomp[2] * jb;
+
+ TRACE(("Raw k1 = %f, k2 = %f, k3 = %f, raw ss = %f\n",k1, k2, k3, (k1 - k3)/k2))
+
+#ifdef ENABLE_DDL
+
+ /* Limit ratio of k3 to k1 to stop zero or -ve ss */
+ if (k3 > (k1 * s->ddulimit)) {
+ k3 = k1 * s->ddulimit;
+ TRACE(("k3 limited to %f due to k3:k1 ratio, ss = %f\n",k3,(k1 - k3)/k2))
+ }
+
+ /* See if there is going to be a problem in fwd */
+ if (k3 < -k1 * s->ddllimit/(1.0 - s->ddllimit)) {
+ /* Adjust ss to allow for fwd limitd computation */
+ k3 = -k1 * s->ddllimit/(1.0 - s->ddllimit);
+ TRACE(("k3 set to %f to allow for fk3:fk2 fwd limit\n",k3))
+ }
+
+#endif /* ENABLE_DDL */
+
+#ifdef DISABLE_TTD
+
+ ss = k1/k2;
+
+#else /* !DISABLE_TTD */
+
+ /* Compute the ab scale factor */
+ ss = (k1 - k3)/k2;
+
+#endif /* !DISABLE_TTD */
+
+#ifdef ENABLE_SS
+ if (ss < SSLLIMIT)
+ ss = SSLLIMIT;
+ else if (ss > SSULIMIT)
+ ss = SSULIMIT;
+#endif /* ENABLE_SS */
+
+ /* Unscale a and b */
+ a = ja / ss;
+ b = jb / ss;
+
+ TRACE(("ss = %f, ttA = %f, a = %f, b = %f\n",ss,ttA,a,b))
+
+ /* Solve for post-adapted cone response of sample */
+ /* using inverse matrix on ttA, a, b */
+ rgba[0] = (20.0/61.0) * ttA
+ + ((41.0 * 11.0)/(61.0 * 23.0)) * a
+ + ((288.0 * 1.0)/(61.0 * 23.0)) * b;
+ rgba[1] = (20.0/61.0) * ttA
+ - ((81.0 * 11.0)/(61.0 * 23.0)) * a
+ - ((261.0 * 1.0)/(61.0 * 23.0)) * b;
+ rgba[2] = (20.0/61.0) * ttA
+ - ((20.0 * 11.0)/(61.0 * 23.0)) * a
+ - ((20.0 * 315.0)/(61.0 * 23.0)) * b;
+
+ TRACE(("rgba %f %f %f\n",rgba[0], rgba[1], rgba[2]))
+
+#ifdef DISABLE_MATRIX
+
+ XYZ[0] = rgba[0];
+ XYZ[1] = rgba[1];
+ XYZ[2] = rgba[2];
+
+#else /* !DISABLE_MATRIX */
+
+#ifdef DISABLE_NONLIN
+ rgbp[0] = 27.13/400.0 * rgba[0];
+ rgbp[1] = 27.13/400.0 * rgba[1];
+ rgbp[2] = 27.13/400.0 * rgba[2];
+#else /* !DISABLE_NONLIN */
+
+ /* Hunt-Pointer_Estevez cone space */
+ /* (with linear segment at the +ve end) */
+ for (i = 0; i < 3; i++) {
+ if (rgba[i] < s->nldxval) {
+ rgbp[i] = s->nldlimit + (rgba[i] - s->nldxval)/s->nldxslope;
+ } else if (rgba[i] <= s->nluxval) {
+ tt = rgba[i] - 0.1;
+ rgbp[i] = pow((27.13 * tt)/(400.0 - tt), 1.0/0.42)/s->Fl;
+// rgbp[i] = pow((27.13 * tt)/(400.0 - tt), 1.0/1.0)/s->Fl;
+ } else {
+ rgbp[i] = s->nlulimit + (rgba[i] - s->nluxval)/s->nluxslope;
+ }
+ }
+#endif /* !DISABLE_NONLIN */
+
+ TRACE(("rgbp %f %f %f\n",rgbp[0], rgbp[1], rgbp[2]))
+
+#ifdef ENABLE_BLUE_ANGLE_FIX
+ ss = rgbp[0] + rgbp[1] + rgbp[2];
+ if (ss < 1e-9)
+ ss = 0.0;
+ else {
+ ss = (rgbp[2]/ss - 1.0/3.0) * 3.0/2.0;
+ if (ss > 0.0)
+ ss = BLUE_BL_MAX * pow(ss, BLUE_BL_POW);
+ }
+ if (ss < 0.0)
+ ss = 0.0;
+ else if (ss > 1.0)
+ ss = 1.0;
+ tt = 0.5 * (rgbp[0] + rgbp[1]);
+ rgbp[0] = (rgbp[0] - ss * tt)/(1.0 - ss);
+ rgbp[1] = (rgbp[1] - ss * tt)/(1.0 - ss);
+
+ TRACE(("rgbp after blue fix %f = %f %f %f\n",ss,rgbp[0], rgbp[1], rgbp[2]))
+#endif
+
+
+#ifdef ENABLE_COMPR
+ /* Undo soft limiting */
+ {
+ double tt; /* Temporary */
+ double wrgb[3]; /* White target */
+
+ /* Make white target white point with same Y value */
+ tt = rgbp[0] * s->icc[1][0] + rgbp[1] * s->icc[1][1] + rgbp[2] * s->icc[1][2];
+ tt = tt > BC_WHMINY ? tt : BC_WHMINY; /* Limit to minimum Y */
+ icmScale3(wrgb, s->rgbpW, tt/s->Wxyz[1]); /* White target at same Y */
+ TRACE(("wrgb %f %f %f\n", wrgb[0], wrgb[1], wrgb[2]))
+
+ /* Un-limit b,g,r in turn */
+ for (i = 2; i >= 0; i--) {
+ double cv; /* Compression value */
+ double ctv; /* Compression target value */
+ double cd; /* Compression displacement needed */
+ double cvec[3]; /* Normalized correction vector */
+ double isec[3]; /* Intersection with plane */
+ double offs; /* Offset of point from orgin*/
+ double range; /* Threshold to start compression */
+ double asym; /* Compression asymtope */
+
+ /* Compute compression direction as vector towards white */
+ /* (We did try correcting in a blend of limit plane normal and white, */
+ /* but compressing towards white seems to be the best.) */
+ icmSub3(cvec, wrgb, rgbp); /* Direction of white target */
+
+ TRACE(("rgbp %f %f %f\n", rgbp[0], rgbp[1], rgbp[2]))
+ TRACE(("cvec %f %f %f\n", cvec[0], cvec[1], cvec[2]))
+
+ if (cvec[i] < 1e-9) { /* compression direction can't correct this coord */
+ TRACE(("Intersection with limit plane failed\n"))
+ continue;
+ }
+
+ /* Scale compression vector to make it move a unit in normal direction */
+ icmScale3(cvec, cvec, 1.0/cvec[i]); /* Normalized vector to white */
+ TRACE(("cvec %f %f %f\n", cvec[0], cvec[1], cvec[2]))
+
+ /* Compute intersection of correction direction with this limit plane */
+ /* (This corresponds with finding displacement of rgbp by cvec */
+ /* such that the current coord value = 0) */
+ icmScale3(isec, cvec, -rgbp[i]); /* (since cvec[i] == 1.0) */
+ icmAdd3(isec, isec, rgbp);
+ TRACE(("isec %f %f %f\n", isec[0], isec[1], isec[2]))
+
+ /* Compute distance from intersection to origin */
+ offs = pow(icmNorm3(isec), 0.85);
+
+ range = s->crange[i] * offs; /* Scale range by distance to origin */
+ if (range > BC_MAXRANGE) /* so that it tapers down as we approach it */
+ range = BC_MAXRANGE; /* and limit maximum */
+
+ /* Aiming above plane when far from origin, */
+ /* but below plane at the origin, so that black isn't affected. */
+ asym = range - 0.2 * (range + (0.01 * s->crange[i]));
+
+ ctv = cv = rgbp[i]; /* Distance above/below limit plane */
+
+ TRACE(("ch %d, offs %f, range %f asym %f, cv %f\n",i, offs,range,asym,cv))
+
+ if (ctv < (range - 1e-12)) { /* Need to expand */
+
+ if (ctv <= asym) {
+ cd = BC_LIMIT;
+ TRACE(("ctv %f < asym %f\n",ctv,asym))
+ } else {
+ double aa, bb;
+ aa = 1.0/(range - ctv);
+ bb = 1.0/(range - asym);
+ if (aa > (bb + 1e-12))
+ cv = range - 1.0/(aa - bb);
+ cd = ctv - cv; /* Displacement needed */
+ }
+ if (cd > BC_LIMIT)
+ cd = BC_LIMIT;
+ TRACE(("ch %d cd = %f, scaled cd %f\n",i,cd,cd))
+
+ if (cd > 1e-9) {
+ icmScale3(cvec, cvec, -cd); /* Compression vector */
+ icmAdd3(rgbp, rgbp, cvec); /* Compress by displacement */
+ TRACE(("rgbp after decomp. = %f %f %f\n",rgbp[0], rgbp[1], rgbp[2]))
+ }
+ }
+ }
+ }
+#endif /* ENABLE_COMPR */
+
+ /* Chromaticaly transformed sample value */
+ /* Spectrally sharpened cone responses */
+ /* XYZ values */
+ icmMulBy3x3(xyz, s->icc, rgbp);
+
+ TRACE(("XYZ = %f %f %f\n",xyz[0], xyz[1], xyz[2]))
+
+ /* Subtract flare */
+ XYZ[0] = s->Fisc * (xyz[0] - s->Fsxyz[0]);
+ XYZ[1] = s->Fisc * (xyz[1] - s->Fsxyz[1]);
+ XYZ[2] = s->Fisc * (xyz[2] - s->Fsxyz[2]);
+
+#endif /* !DISABLE_MATRIX */
+
+ TRACE(("XYZ after flare = %f %f %f\n",XYZ[0], XYZ[1], XYZ[2]))
+ TRACE(("\n"))
+
+#ifdef DIAG2
+ printf("Processing:\n");
+ printf("Jab = %f %f %f\n", Jabi[0], Jabi[1], Jabi[2]);
+ printf("Chroma C = %f\n", C);
+ printf("Preliminary Saturation ss = %f\n", ss);
+ printf("Lightness J = %f, H.K. Lightness = %f\n", J * 100, JJ * 100);
+ printf("Achromatic response A = %f\n", A);
+ printf("Eccentricity factor e = %f\n", e);
+ printf("Hue angle h = %f\n", h);
+ printf("Post adapted cone response rgba = %f %f %f\n", rgba[0], rgba[1], rgba[2]);
+ printf("Hunundeft-P-E cone space rgbp = %f %f %f\n", rgbp[0], rgbp[1], rgbp[2]);
+ printf("Including flare XYZ = %f %f %f\n", xyz[0], xyz[1], xyz[2]);
+ printf("XYZ = %f %f %f\n", XYZ[0], XYZ[1], XYZ[2]);
+ printf("\n");
+#endif
+ return 0;
+}
+