Initial import of argyll version 1.5.1-8debian/1.5.1-8

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;
+}
+