/* 
 * International Color Consortium color transform expanded support
 *
 * Author:  Graeme W. Gill
 * Date:    2/7/00
 * Version: 1.00
 *
 * Copyright 2000, 2001 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.
 *
 * Based on the old iccXfm class.
 */

/*
 * This module expands the basic icclib functionality,
 * providing more functionality in exercising.
 * The implementation for the three different types
 * of profile representation, are in their own source files.
 */

/*
 * TTBD:
 *       Some of the error handling is crude. Shouldn't use
 *       error(), should return status.
 *
 */

#include <sys/types.h>
#include <string.h>
#include <ctype.h>
#ifdef __sun
#include <unistd.h>
#endif
#if defined(__IBMC__) && defined(_M_IX86)
#include <float.h>
#endif
#include "numlib.h"
#include "counters.h"
#include "plot.h"
#include "../h/sort.h"
#include "xicc.h"		/* definitions for this library */

#define USE_CAM			/* Use CIECAM02 for clipping and gamut mapping, else use Lab */

#undef DEBUG			/* Plot 1d Luts */

#ifdef DEBUG
#include "plot.h"
#endif

#define MAX_INVSOLN 4

static void xicc_del(xicc *p);
icxLuBase * xicc_get_luobj(xicc *p, int flags, icmLookupFunc func, icRenderingIntent intent,
                           icColorSpaceSignature pcsor, icmLookupOrder order,
                           icxViewCond *vc, icxInk *ink);
static icxLuBase *xicc_set_luobj(xicc *p, icmLookupFunc func, icRenderingIntent intent,
                            icmLookupOrder order, int flags, int no, int nobw, cow *points,
							icxMatrixModel *skm,
                            double dispLuminance, double wpscale, double smooth, double avgdev,
                            icxViewCond *vc, icxInk *ink, xcal *cal, int quality);
static void icxLutSpaces(icxLuBase *p, icColorSpaceSignature *ins, int *inn,
                         icColorSpaceSignature *outs, int *outn,
                         icColorSpaceSignature *pcs);
static void icxLuSpaces(icxLuBase *p, icColorSpaceSignature *ins, int *inn,
                        icColorSpaceSignature *outs, int *outn, 
                        icmLuAlgType *alg, icRenderingIntent *intt,
                        icmLookupFunc *fnc, icColorSpaceSignature *pcs);
static void icxLu_get_native_ranges (icxLuBase *p,
                              double *inmin, double *inmax, double *outmin, double *outmax);
static void icxLu_get_ranges (icxLuBase *p,
                              double *inmin, double *inmax, double *outmin, double *outmax);
static void icxLuEfv_wh_bk_points(icxLuBase *p, double *wht, double *blk, double *kblk);
int xicc_get_viewcond(xicc *p, icxViewCond *vc);

/* The different profile types are in their own source filesm */
/* and are included to keep their functions private. (static) */
#include "xmono.c"
#include "xmatrix.c"
#include "xlut.c"		/* New xfit3 in & out optimising based profiles */
//#include "xlut1.c"	/* Old xfit1 device curve based profiles */

#ifdef NT		/* You'd think there might be some standards.... */
# ifndef __BORLANDC__
#  define stricmp _stricmp
# endif
#else
# define stricmp strcasecmp
#endif
/* Utilities */

/* Return a string description of the given enumeration value */
const char *icx2str(icmEnumType etype, int enumval) {

	if (etype == icmColorSpaceSignature) {
		if (((icColorSpaceSignature)enumval) == icxSigJabData)
			return "Jab";
		else if (((icColorSpaceSignature)enumval) == icxSigJChData)
			return "JCh";
		else if (((icColorSpaceSignature)enumval) == icxSigLChData)
			return "LCh";
	} else if (etype == icmRenderingIntent) {
		if (((icRenderingIntent)enumval) == icxAppearance)
			return "icxAppearance";
		else if (((icRenderingIntent)enumval) == icxAbsAppearance)
			return "icxAbsAppearance";
		else if (((icRenderingIntent)enumval) == icxPerceptualAppearance)
			return "icxPerceptualAppearance";
		else if (((icRenderingIntent)enumval) == icxAbsPerceptualAppearance)
			return "icxAbsPerceptualAppearance";
		else if (((icRenderingIntent)enumval) == icxSaturationAppearance)
			return "icxSaturationAppearance";
		else if (((icRenderingIntent)enumval) == icxAbsSaturationAppearance)
			return "icxAbsSaturationAppearance";
	}
	return icm2str(etype, enumval);
}

/* Common xicc stuff */

/* Return information about the native lut in/out colorspaces. */
/* Any pointer may be NULL if value is not to be returned */
static void
icxLutSpaces(
	icxLuBase *p,					/* This */
	icColorSpaceSignature *ins,		/* Return input color space */
	int *inn,						/* Return number of input components */
	icColorSpaceSignature *outs,	/* Return output color space */
	int *outn,						/* Return number of output components */
	icColorSpaceSignature *pcs		/* Return PCS color space */
) {
	p->plu->lutspaces(p->plu, ins, inn, outs, outn, pcs);
}

/* Return information about the overall lookup in/out colorspaces, */
/* including allowance for any PCS override. */
/* Any pointer may be NULL if value is not to be returned */
static void
icxLuSpaces(
	icxLuBase *p,					/* This */
	icColorSpaceSignature *ins,		/* Return input color space */
	int *inn,						/* Return number of input components */
	icColorSpaceSignature *outs,	/* Return output color space */
	int *outn,						/* Return number of output components */
	icmLuAlgType *alg,				/* Return type of lookup algorithm used */
    icRenderingIntent *intt,		/* Return the intent implemented */
    icmLookupFunc *fnc,				/* Return the profile function being implemented */
	icColorSpaceSignature *pcs		/* Return the effective PCS */
) {
    icmLookupFunc function;
	icColorSpaceSignature npcs;		/* Native PCS */

	p->plu->spaces(p->plu, NULL, inn, NULL, outn, alg, NULL, &function, &npcs, NULL);

	if (intt != NULL)
		*intt = p->intent;

	if (fnc != NULL)
		*fnc = function;

	if (ins != NULL)
		*ins = p->ins;

	if (outs != NULL)
		*outs = p->outs;

	if (pcs != NULL)
		*pcs = p->pcs;
}

/* Return the native (internaly visible) colorspace value ranges */
static void
icxLu_get_native_ranges (
icxLuBase *p,
double *inmin, double *inmax,		/* Return maximum range of inspace values */
double *outmin, double *outmax		/* Return maximum range of outspace values */
) {
	int i;
	if (inmin != NULL) {
		for (i = 0; i < p->inputChan; i++)
			inmin[i] = p->ninmin[i];
	}
	if (inmax != NULL) {
		for (i = 0; i < p->inputChan; i++)
			inmax[i] = p->ninmax[i];
	}
	if (outmin != NULL) {
		for (i = 0; i < p->outputChan; i++)
			outmin[i] = p->noutmin[i];
	}
	if (outmax != NULL) {
		for (i = 0; i < p->outputChan; i++)
			outmax[i] = p->noutmax[i];
	}
}

/* Return the effective (externaly visible) colorspace value ranges */
static void
icxLu_get_ranges (
icxLuBase *p,
double *inmin, double *inmax,		/* Return maximum range of inspace values */
double *outmin, double *outmax		/* Return maximum range of outspace values */
) {
	int i;
	if (inmin != NULL) {
		for (i = 0; i < p->inputChan; i++)
			inmin[i] = p->inmin[i];
	}
	if (inmax != NULL) {
		for (i = 0; i < p->inputChan; i++)
			inmax[i] = p->inmax[i];
	}
	if (outmin != NULL) {
		for (i = 0; i < p->outputChan; i++)
			outmin[i] = p->outmin[i];
	}
	if (outmax != NULL) {
		for (i = 0; i < p->outputChan; i++)
			outmax[i] = p->outmax[i];
	}
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Routine to figure out a suitable black point for CMYK */

/* Structure to hold optimisation information */
typedef struct {
	icmLuBase *p;
	int kch;				/* K channel, -1 if none */
	double tlimit, klimit;	/* Ink limit values */
	int inn;				/* Number of input channels */
	icColorSpaceSignature outs;		/* Output space */
	double p1[3];			/* white pivot point in abs Lab */
	double p2[3];			/* Point on vector towards black in abs Lab */
	double toll;			/* Tollerance of black direction */
} bpfind;

/* Optimise device values to minimise L, while remaining */
/* within the ink limit, and staying in line between p1 (white) and p2 (black dir) */
static double bpfindfunc(void *adata, double pv[]) {
	bpfind *b = (bpfind *)adata;
	double rv = 0.0;
	double Lab[3];
	double lr, ta, tb, terr;	/* L ratio, target a, target b, target error */
	double ovr = 0.0;
	int e;

	/* Compute amount outside total limit */
	if (b->tlimit >= 0.0) {
		double sum;
		for (sum = 0.0, e = 0; e < b->inn; e++)
			sum += pv[e];
		if (sum > b->tlimit) {
			ovr = sum - b->tlimit;
#ifdef DEBUG
	printf("~1 total ink ovr = %f\n",ovr);
#endif
		}
	}

	/* Compute amount outside black limit */
	if (b->klimit >= 0.0 && b->kch >= 0) {
		double kval = pv[b->kch] - b->klimit;
		if (kval > ovr) {
			ovr = kval;
#ifdef DEBUG
	printf("~1 black ink ovr = %f\n",ovr);
#endif
		}
	}
	/* Compute amount outside device value limits 0.0 - 1.0 */
	{
		double dval;
		for (dval = -1.0, e = 0; e < b->inn; e++) {
			if (pv[e] < 0.0) {
				if (-pv[e] > dval)
					dval = -pv[e];
			} else if (pv[e] > 1.0) {
				if ((pv[e] - 1.0) > dval)
					dval = pv[e] - 1.0;
			}
		}
		if (dval > ovr)
			ovr = dval;
	}

	/* Compute the Lab value: */
	b->p->lookup(b->p, Lab, pv);
	if (b->outs == icSigXYZData)
		icmXYZ2Lab(&icmD50, Lab, Lab);

#ifdef DEBUG
	printf("~1 p1 =  %f %f %f, p2 = %f %f %f\n",b->p1[0],b->p1[1],b->p1[2],b->p2[0],b->p2[1],b->p2[2]);
	printf("~1 device value %f %f %f %f, Lab = %f %f %f\n",pv[0],pv[1],pv[2],pv[3],Lab[0],Lab[1],Lab[2]);
#endif

	/* Primary aim is to minimise L value */
	rv = Lab[0];

	/* See how out of line from p1 to p2 we are */
	lr = (Lab[0] - b->p1[0])/(b->p2[0] - b->p1[0]);		/* Distance towards p2 from p1 */
	ta = lr * (b->p2[1] - b->p1[1]) + b->p1[1];			/* Target a value */
	tb = lr * (b->p2[2] - b->p1[2]) + b->p1[2];			/* Target b value */

	terr = (ta - Lab[1]) * (ta - Lab[1])
	     + (tb - Lab[2]) * (tb - Lab[2]);
	
	if (terr < b->toll)		/* Tollerance error doesn't count until it's over tollerance */
		terr = 0.0;
	
#ifdef DEBUG
	printf("~1 target error %f\n",terr);
#endif
	rv += XICC_BLACK_FIND_ABERR_WEIGHT * terr;	/* Make ab match more important than min. L */

#ifdef DEBUG
	printf("~1 out of range error %f\n",ovr);
#endif
	rv += 200 * ovr;

#ifdef DEBUG
	printf("~1 black find tc ret %f\n",rv);
#endif
	return rv;
}

/* Try and compute a real black point in XYZ given an iccLu, */
/* and also return the K only black or the normal black if the device doesn't have K */
/* black[] will be unchanged if black cannot be computed. */
/* Note that the black point will be in the space of the Lu */
/* converted to XYZ, so will have the Lu's intent etc. */
/* (Note that this is duplicated in xlut.c set_icxLuLut() !!!) */
static void icxLu_comp_bk_point(
icxLuBase *x,
int gblk,				/* If nz, compute black if possible. */
double *white,			/* XYZ Input, used for computing black */
double *black,			/* XYZ Input & Output. Set if gblk NZ and can be computed */
double *kblack			/* XYZ Output. Looked up if possible or set to black[] otherwise */
) {
	icmLuBase *p = x->plu;
	icmLuBase *op = p;			/* Original icmLu, in case we replace p */
	icc *icco = p->icp;
	icmHeader *h = icco->header;
	icColorSpaceSignature ins, outs;
	int inn, outn;
	icmLuAlgType alg;
	icRenderingIntent intt;
	icmLookupFunc fnc;
	icmLookupOrder ord;
	int kch = -1;
	double dblack[MAX_CHAN];	/* device black value */
	int e;

#ifdef DEBUG
	printf("~1 icxLu_comp_bk_point() called, gblk %d, white = %s, black = %s\n",gblk,icmPdv(3, white),icmPdv(3,black));
#endif
	/* Default return incoming black as K only black */
	kblack[0] = black[0];
	kblack[1] = black[1];
	kblack[2] = black[2];

	/* Get the effective characteristics of the Lu */
	p->spaces(p, &ins, &inn, &outs, &outn, &alg, &intt, &fnc, NULL, &ord);

	if (fnc == icmBwd) { /* Hmm. We've got PCS to device, and we want device to PCS. */

		/* Strictly speaking this is a dubious approach, since for a cLut profile */
		/* the B2A table could make the effective white and black points */
		/* anything it likes, and they don't have to match what the corresponding */
		/* A2B table does. In our usage it's probably OK, since we tend */
		/* to use colorimetric B2A */ 
#ifdef DEBUG
		printf("~1 getting icmFwd\n");
#endif
		if ((p = icco->get_luobj(icco, icmFwd, intt, ins, ord)) == NULL)
			error("icxLu_comp_bk_point: assert: getting Fwd Lookup failed!");

		p->spaces(p, &ins, &inn, &outs, &outn, &alg, &intt, &fnc, NULL, &ord);
	}

	if (outs != icSigXYZData && outs != icSigLabData) {
		error("icxLu_comp_bk_point: assert: icc Lu output is not XYZ or Lab!, outs = 0x%x, ");
	}

#ifdef DEBUG
	printf("~1 icxLu_comp_bk_point called for inn = %d, ins = %s\n", inn, icx2str(icmColorSpaceSignature,ins));
#endif

	switch (ins) {

    	case icSigXYZData:
    	case icSigLabData:
    	case icSigLuvData:
    	case icSigYxyData:
#ifdef DEBUG
			printf("~1 Assuming CIE colorspace black is 0.0\n");
#endif
			if (gblk) {
				for (e = 0; e < inn; e++)
					black[0] = 0.0;
			}
			kblack[0] = black[0];
			kblack[1] = black[1];
			kblack[2] = black[2];
			return;

		case icSigRgbData:
#ifdef DEBUG
			printf("~1 RGB:\n");
#endif
			for (e = 0; e < inn; e++)
				dblack[e] = 0.0;
			break;

		case icSigGrayData: {	/* Could be additive or subtractive */
			double dval[1];
			double minv[3], maxv[3];
#ifdef DEBUG
			printf("~1 Gray:\n");
#endif
			/* Check out 0 and 100% colorant */
			dval[0] = 0.0;
			p->lookup(p, minv, dval);
			if (outs == icSigXYZData)
				icmXYZ2Lab(&icmD50, minv, minv);
			dval[0] = 1.0;
			p->lookup(p, maxv, dval);
			if (outs == icSigXYZData)
				icmXYZ2Lab(&icmD50, maxv, maxv);

			if (minv[0] < maxv[0])
				dblack[0] = 0.0;
			else
				dblack[0] = 1.0;
			}
			break;

		case icSigCmyData:
			for (e = 0; e < inn; e++)
				dblack[e] = 1.0;
			break;

		case icSigCmykData:
#ifdef DEBUG
			printf("~1 CMYK:\n");
#endif
			kch = 3;
			dblack[0] = 0.0;
			dblack[1] = 0.0;
			dblack[2] = 0.0;
			dblack[3] = 1.0;
			if (alg == icmLutType) {
				icxLuLut *pp = (icxLuLut *)x;
		
				if (pp->ink.tlimit >= 0.0)
					dblack[kch] = pp->ink.tlimit;
			};
			break;

			/* Use a heursistic. */
			/* This duplicates code in icxGetLimits() :-( */
			/* Colorant guessing should go in icclib ? */
		case icSig2colorData:
    	case icSig3colorData:
    	case icSig4colorData:
    	case icSig5colorData:
    	case icSig6colorData:
    	case icSig7colorData:
    	case icSig8colorData:
    	case icSig9colorData:
    	case icSig10colorData:
    	case icSig11colorData:
    	case icSig12colorData:
    	case icSig13colorData:
    	case icSig14colorData:
    	case icSig15colorData:
    	case icSigMch5Data:
    	case icSigMch6Data:
    	case icSigMch7Data:
    	case icSigMch8Data: {
				double dval[MAX_CHAN];
				double ncval[3];
				double cvals[MAX_CHAN][3];
				int nlighter, ndarker;

				/* Decide if the colorspace is additive or subtractive */
#ifdef DEBUG
				printf("~1 N channel:\n");
#endif

				/* First the no colorant value */
				for (e = 0; e < inn; e++)
					dval[e] = 0.0;
				p->lookup(p, ncval, dval);
				if (outs == icSigXYZData)
					icmXYZ2Lab(&icmD50, ncval, ncval);

				/* Then all the colorants */
				nlighter = ndarker = 0;
				for (e = 0; e < inn; e++) {
					dval[e] = 1.0;
					p->lookup(p, cvals[e], dval);
					if (outs == icSigXYZData)
						icmXYZ2Lab(&icmD50, cvals[e], cvals[e]);
					dval[e] = 0.0;
					if (fabs(cvals[e][0] - ncval[0]) > 5.0) {
						if (cvals[e][0] > ncval[0])
							nlighter++;
						else
							ndarker++;
					}
				}
				if (ndarker == 0 && nlighter > 0) {		/* Assume additive */
					for (e = 0; e < inn; e++)
						dblack[e] = 0.0;
#ifdef DEBUG
					printf("~1 N channel is additive:\n");
#endif

				} else if (ndarker > 0 && nlighter == 0) {				/* Assume subtractive. */
					double pbk[3] = { 0.0,0.0,0.0 };	/* Perfect black */
					double smd = 1e10;			/* Smallest distance */

#ifdef DEBUG
					printf("~1 N channel is subtractive:\n");
#endif
					/* See if we can guess the black channel */
					for (e = 0; e < inn; e++) {
						double tt;
						tt = icmNorm33sq(pbk, cvals[e]);
						if (tt < smd) {
							smd = tt;	
							kch = e;
						}
					}
					/* See if the black seems sane */
					if (cvals[kch][0] > 40.0
					 || fabs(cvals[kch][1]) > 10.0
					 || fabs(cvals[kch][2]) > 10.0) {
						if (p != op)
							p->del(p);
#ifdef DEBUG
						printf("~1 black doesn't look sanem so assume nothing\n");
#endif
						return;		/* Assume nothing */
					}

					/* Chosen kch as black */
					for (e = 0; e < inn; e++)
						dblack[e] = 0.0;
					dblack[kch] = 1.0;
					if (alg == icmLutType) {
						icxLuLut *pp = (icxLuLut *)x;
			
						if (pp->ink.tlimit >= 0.0)
							dblack[kch] = pp->ink.tlimit;
					};
#ifdef DEBUG
					printf("~1 N channel K = chan %d\n",kch);
#endif
				} else {
					if (p != op)
						p->del(p);
#ifdef DEBUG
					printf("~1 can't figure if additive or subtractive, so assume nothing\n");
#endif
					return;			/* Assume nothing */
				}
			}
			break;

		default:
#ifdef DEBUG
			printf("~1 unhandled colorspace, so assume nothing\n");
#endif
			if (p != op)
				p->del(p);
			return;				/* Don't do anything */
	}

	/* Lookup the K only value */
	if (kch >= 0) {
		p->lookup(p, kblack, dblack);
		
		/* We always return XYZ */
		if (outs == icSigLabData)
			icmLab2XYZ(&icmD50, kblack, kblack);
	}

	if (gblk == 0) {		/* That's all we have to do */
#ifdef DEBUG
		printf("~1 gblk == 0, so only return kblack\n");
#endif
		if (p != op)
			p->del(p);
		return;
	}

	/* Lookup the device black or K only value as a default */
	p->lookup(p, black, dblack);		/* May be XYZ or Lab */
#ifdef DEBUG
	printf("~1 Got default lu black %f %f %f, kch = %d\n", black[0],black[1],black[2],kch);
#endif

	/* !!! Hmm. For CMY and RGB we are simply using the device */
	/* combination values as the black point. In reality we might */
	/* want to have the option of using a neutral black point, */
	/* just like CMYK ?? */

	if (kch >= 0) {	/* The space is subtractive with a K channel. */
		/* If XICC_NEUTRAL_CMYK_BLACK then locate the darkest */
		/* CMYK within limits with the same chromaticity as the white point, */
		/* otherwise locate the device value within the ink limits that is */
		/* in the direction of the K channel */
		bpfind bfs;					/* Callback context */
		double sr[MXDO];			/* search radius */
		double tt[MXDO];			/* Temporary */
		double rs0[MXDO], rs1[MXDO];	/* Random start candidates */
		int trial;
		double brv;

		/* Setup callback function context */
		bfs.p = p;
		bfs.inn = inn;
		bfs.outs = outs;

		bfs.kch = kch;
		bfs.tlimit = -1.0;
		bfs.klimit = -1.0;
		bfs.toll = XICC_BLACK_POINT_TOLL;

		if (alg == icmLutType) {
			icxLuLut *pp = (icxLuLut *)x;

			pp->kch = kch;
			bfs.tlimit = pp->ink.tlimit;
			bfs.klimit = pp->ink.klimit;
#ifdef DEBUG
			printf("~1 tlimit = %f, klimit = %f\n",bfs.tlimit,bfs.klimit);
#endif
		};
	
#ifdef XICC_NEUTRAL_CMYK_BLACK
#ifdef DEBUG
		printf("~1 Searching for neutral black\n");
#endif
		/* white has been given to us in XYZ */
		icmXYZ2Lab(&icmD50, bfs.p1, white);		/* pivot Lab */
		icmCpy3(bfs.p2, white);					/* temp white XYZ */
		icmScale3(bfs.p2, bfs.p2, 0.02);	 /* Scale white XYZ towards 0,0,0 */
		icmXYZ2Lab(&icmD50, bfs.p2, bfs.p2); /* Convert black direction to Lab */

#else /* Use K directin black */
#ifdef DEBUG
		printf("~1 Searching for K direction black\n");
#endif
		icmXYZ2Lab(&icmD50, bfs.p1, white);			/* Pivot */

		/* Now figure abs Lab value of K only, as the direction */
		/* to use for the rich black. */
		for (e = 0; e < inn; e++)
			dblack[e] = 0.0;
		if (bfs.klimit < 0.0)
			dblack[kch] = 1.0;
		else
			dblack[kch] = bfs.klimit;		/* K value */

		p->lookup(p, black, dblack);

		if (outs == icSigXYZData) {
			icmXYZ2Lab(&icmD50, bfs.p2, black);		/* K direction */
		} else {
			icmAry2Ary(bfs.p2, black);
		}
#endif

#ifdef DEBUG
		printf("~1 Lab pivot %f %f %f, Lab K direction %f %f %f\n",bfs.p1[0],bfs.p1[1],bfs.p1[2],bfs.p2[0],bfs.p2[1],bfs.p2[2]);
#endif
		/* Start with the K only as the current best value */
		brv = bpfindfunc((void *)&bfs, dblack);
#ifdef DEBUG
		printf("~1 initial brv for K only = %f\n",brv);
#endif

		/* Set the random start 0 location as 000K */
		/* and the random start 1 location as CMY0 */
		{
			double tt;

			for (e = 0; e < inn; e++)
				dblack[e] = rs0[e] = 0.0;
			if (bfs.klimit < 0.0)
				dblack[kch] = rs0[kch] = 1.0;
			else
				dblack[kch] = rs0[kch] = bfs.klimit;		/* K value */

			if (bfs.tlimit < 0.0)
				tt = 1.0;
			else
				tt = bfs.tlimit/(inn - 1.0);
			for (e = 0; e < inn; e++)
				rs1[e] = tt;
			rs1[kch] = 0.0;		/* K value */
		}

		/* Find the device black point using optimization */
		/* Do several trials to avoid local minima. */
		rand32(0x12345678);	/* Make trial values deterministic */
		for (trial = 0; trial < 200; trial++) {
			double rv;			/* Temporary */

			/* Start first trial at 000K */
			if (trial == 0) {
				for (e = 0; e < inn; e++) {
					tt[e] = rs0[e];
					sr[e] = 0.1;
				}

			} else {
				/* Base is random between 000K and CMY0: */
				if (trial < 100) {
					rv = d_rand(0.0, 1.0);
					for (e = 0; e < inn; e++) {
						tt[e] = rv * rs0[e] + (1.0 - rv) * rs1[e];
						sr[e] = 0.1;
					}
				/* Base on current best */
				} else {
					for (e = 0; e < inn; e++) {
						tt[e] = dblack[e];
						sr[e] = 0.1;
					}
				}
	
				/* Then add random start offset */
				for (rv = 0.0, e = 0; e < inn; e++) {
					tt[e] +=  d_rand(-0.5, 0.5);
					if (tt[e] < 0.0)
						tt[e] = 0.0;
					else if (tt[e] > 1.0)
						tt[e] = 1.0;
				}
			}

			/* Clip black */
			if (bfs.klimit >= 0.0 && tt[kch] > bfs.klimit)
				tt[kch] = bfs.klimit;	
			
			/* Compute amount outside total limit */
			if (bfs.tlimit >= 0.0) {
				for (rv = 0.0, e = 0; e < inn; e++)
					rv += tt[e];
		
				if (rv > bfs.tlimit) {
					rv /= (double)inn;
					for (e = 0; e < inn; e++)
						tt[e] -= rv;
				}
			}
	
			if (powell(&rv, inn, tt, sr, 0.000001, 1000, bpfindfunc,
			                      (void *)&bfs, NULL, NULL) == 0) {
#ifdef DEBUG
				printf("~1 trial %d, rv %f bp %f %f %f %f\n",trial,rv,tt[0],tt[1],tt[2],tt[3]);
#endif
				if (rv < brv) {
#ifdef DEBUG
					printf("~1   new best\n");
#endif
					brv = rv;
					for (e = 0; e < inn; e++)
						dblack[e] = tt[e];
				}
			}
		}
		if (brv > 1000.0)
			error("icxLu_comp_bk_point: Black point powell failed");

		for (e = 0; e < inn; e++) { /* Make sure device values are in range */
			if (dblack[e] < 0.0)
				dblack[e] = 0.0;
			else if (dblack[e] > 1.0)
				dblack[e] = 1.0;
		}
		/* Now have device black in dblack[] */
#ifdef DEBUG
		printf("~1 got device black %f %f %f %f\n",dblack[0], dblack[1], dblack[2], dblack[3]);
#endif

		p->lookup(p, black, dblack);		/* Convert to PCS */
	}

	if (p != op)
		p->del(p);

	/* We always return XYZ */
	if (outs == icSigLabData)
		icmLab2XYZ(&icmD50, black, black);

#ifdef DEBUG
	printf("~1 returning %f %f %f\n", black[0], black[1], black[2]);
#endif

	return;
}

/* - - - - - - - - - - - - - - - - - - - - - - - */

/* Return the media white and black points */
/* in the xlu effective PCS colorspace. Pointers may be NULL. */
/* (ie. these will be relative values for relative intent etc.) */
static void icxLuEfv_wh_bk_points(
icxLuBase *p,
double *wht,
double *blk,
double *kblk		/* K only black */
) {
	double white[3], black[3], kblack[3];

	/* Get the Lu PCS converted to XYZ icc black and white points in XYZ */
	if (p->plu->lu_wh_bk_points(p->plu, white, black)) {
		/* Black point is assumed. We should determine one instead. */
		/* Lookup K only black too */
		icxLu_comp_bk_point(p, 1, white, black, kblack);

	} else {
		/* Lookup a possible K only black */
		icxLu_comp_bk_point(p, 0, white, black, kblack);
	}

//printf("~1 white %f %f %f, black %f %f %f, kblack %f %f %f\n",white[0],white[1],white[2],black[0],black[1],black[2],kblack[0],kblack[1],kblack[2]);
	/* Convert to possibl xicc override PCS */
	switch (p->pcs) {
		case icSigXYZData:
			break;								/* Don't have to do anyting */
		case icSigLabData:
			icmXYZ2Lab(&icmD50, white, white);	/* Convert from XYZ to Lab */
			icmXYZ2Lab(&icmD50, black, black);
			icmXYZ2Lab(&icmD50, kblack, kblack);
			break;
		case icxSigJabData:
			p->cam->XYZ_to_cam(p->cam, white, white);	/* Convert from XYZ to Jab */
			p->cam->XYZ_to_cam(p->cam, black, black);
			p->cam->XYZ_to_cam(p->cam, kblack, kblack);
			break;
		default:
			break;
	}

//printf("~1 icxLuEfv_wh_bk_points: pcsor %s White %f %f %f, Black %f %f %f\n", icx2str(icmColorSpaceSignature,p->pcs), white[0], white[1], white[2], black[0], black[1], black[2]);
	if (wht != NULL) {
		wht[0] = white[0];
		wht[1] = white[1];
		wht[2] = white[2];
	}

	if (blk != NULL) {
		blk[0] = black[0];
		blk[1] = black[1];
		blk[2] = black[2];
	}

	if (kblk != NULL) {
		kblk[0] = kblack[0];
		kblk[1] = kblack[1];
		kblk[2] = kblack[2];
	}
}

/* Create an instance of an xicc object */
xicc *new_xicc(
icc *picc		/* icc we are expanding */
) {
	xicc *p;
	if ((p = (xicc *) calloc(1,sizeof(xicc))) == NULL)
		return NULL;
	p->pp = picc;
	p->del           = xicc_del;
	p->get_luobj     = xicc_get_luobj;
	p->set_luobj     = xicc_set_luobj;
	p->get_viewcond  = xicc_get_viewcond;

	/* Create an xcal if there is the right tag in the profile */
	p->cal = xiccReadCalTag(p->pp);
	p->nodel_cal = 0;	/* We created it, we will delete it */

	return p;
}

/* Do away with the xicc (but not the icc!) */
static void xicc_del(
xicc *p
) {
	if (p->cal != NULL && p->nodel_cal == 0)
		p->cal->del(p->cal);
	free (p);
}


/* Return an expanded lookup object, initialised */
/* from the icc. */
/* Return NULL on error, check errc+err for reason. */
/* Set the pcsor & intent to consistent and values if */
/* Jab and/or icxAppearance has been requested. */
/* Create the underlying icm lookup object that is used */
/* to create and implement the icx one. The icm will be used */
/* to translate from native to effective PCS, unless the */
/* effective PCS is Jab, in which case the icm will be set to */
/* have an effective PCS of XYZ. Since native<->effecive PCS conversion */
/* is done at the to/from_abs() stage, none of this affects the individual */
/* conversion steps, which will all talk the native PCS (unless merged). */
icxLuBase *xicc_get_luobj(
xicc *p,					/* this */
int flags,					/* clip, merge flags */
icmLookupFunc func,			/* Functionality */
icRenderingIntent intent,	/* Intent */
icColorSpaceSignature pcsor,/* PCS override (0 = def) */
icmLookupOrder order,		/* Search Order */
icxViewCond *vc,			/* Viewing Condition (may be NULL if pcsor is not CIECAM) */
icxInk *ink					/* inking details (NULL for default) */
) {
	icmLuBase *plu;
	icxLuBase *xplu;
	icmLuAlgType alg;
	icRenderingIntent n_intent = intent;			/* Native Intent to request */
	icColorSpaceSignature n_pcs = icmSigDefaultData;	/* Native PCS to request */

//printf("~1 xicc_get_luobj got intent %s and pcsor %s\n",icx2str(icmRenderingIntent,intent),icx2str(icmColorSpaceSignature,pcsor));

	/* Ensure that appropriate PCS is slected for an appearance intent */
	if (intent == icxAppearance
	 || intent == icxAbsAppearance
	 || intent == icxPerceptualAppearance
	 || intent == icxAbsPerceptualAppearance
	 || intent == icxSaturationAppearance
	 || intent == icxAbsSaturationAppearance) {
		pcsor = icxSigJabData;

	/* Translate non-Jab intents to the equivalent appearance "intent" if pcsor == Jab. */
	/* This is how we get these when the UI's don't list all the apperances intents, */
	/* we select the analogous non-apperance intent with pcsor = Jab. */
	/* Note that Abs/non-abs selects between Apperance and AbsAppearance. */
	} else if (pcsor == icxSigJabData) {
		if (intent == icRelativeColorimetric)
			intent = icxAppearance;
		else if (intent == icAbsoluteColorimetric)
			intent = icxAbsAppearance;
		else if (intent == icPerceptual)
			intent = icxPerceptualAppearance;
		else if (intent == icmAbsolutePerceptual)
			intent = icxAbsPerceptualAppearance;
		else if (intent == icSaturation)
			intent = icxSaturationAppearance;
		else if (intent == icmAbsoluteSaturation)
			intent = icxAbsSaturationAppearance;
		else
			intent = icxAppearance;
	}

	/* Translate intent asked for into intent needed in icclib */
	if      (intent == icxAppearance
	      || intent == icxAbsAppearance)
		n_intent = icAbsoluteColorimetric;
	else if (intent == icxPerceptualAppearance
	      || intent == icxAbsPerceptualAppearance)
		n_intent = icmAbsolutePerceptual;
	else if (intent == icxSaturationAppearance
	      || intent == icxAbsSaturationAppearance)
		n_intent = icmAbsoluteSaturation;

	if (pcsor != icmSigDefaultData)
		n_pcs = pcsor;			/* There is an icclib override */

	if (pcsor == icxSigJabData)	/* xicc override */
		n_pcs = icSigXYZData;	/* Translate to XYZ */

//printf("~1 xicc_get_luobj processed intent %s and pcsor %s\n",icx2str(icmRenderingIntent,intent),icx2str(icmColorSpaceSignature,pcsor));
//printf("~1 xicc_get_luobj icclib intent %s and pcsor %s\n",icx2str(icmRenderingIntent,n_intent),icx2str(icmColorSpaceSignature,n_pcs));
	/* Get icclib lookup object */
	if ((plu = p->pp->get_luobj(p->pp, func, n_intent, n_pcs, order)) == NULL) {
		p->errc = p->pp->errc;		/* Copy error */
		strcpy(p->err, p->pp->err);
		return NULL;
	}

	/* Figure out what the algorithm is */
	plu->spaces(plu, NULL, NULL, NULL, NULL, &alg, NULL, NULL, &n_pcs, NULL);

	/* make sure its "Abs CAM" */
	if (vc!= NULL
	 && (intent == icxAbsAppearance
	  || intent == icxAbsPerceptualAppearance
	  || intent == icxAbsSaturationAppearance)) {	/* make sure its "Abs CAM" */
		/* Set white point and flare color to D50 */
		/* (Hmm. This doesn't match what happens within collink with absolute intent!!) */
		vc->Wxyz[0] = icmD50.X/icmD50.Y;
		vc->Wxyz[1] = icmD50.Y/icmD50.Y;	// Normalise white reference to Y = 1 ?
		vc->Wxyz[2] = icmD50.Z/icmD50.Y;
	
		vc->Fxyz[0] = icmD50.X;
		vc->Fxyz[1] = icmD50.Y;
		vc->Fxyz[2] = icmD50.Z;
	}

	/* Call xiccLu wrapper creation */
	switch (alg) {
		case icmMonoFwdType:
			xplu = new_icxLuMono(p, flags, plu, func, intent, pcsor, vc, 0);
			break;
		case icmMonoBwdType:
			xplu = new_icxLuMono(p, flags, plu, func, intent, pcsor, vc, 1);
			break;
    	case icmMatrixFwdType:
			xplu = new_icxLuMatrix(p, flags, plu, func, intent, pcsor, vc, 0);
			break;
    	case icmMatrixBwdType:
			xplu = new_icxLuMatrix(p, flags, plu, func, intent, pcsor, vc, 1);
			break;
    	case icmLutType:
			xplu = new_icxLuLut(p, flags, plu, func, intent, pcsor, vc, ink);
			break;
		default:
			xplu = NULL;
			break;
	}

	return xplu;
}


/* Return an expanded lookup object, initialised */
/* from the icc, and then overwritten by a conversion */
/* created from the supplied scattered data points. */
/* The Lut is assumed to be a device -> native PCS profile. */
/* If the SET_WHITE and/or SET_BLACK flags are set, */
/* discover the white/black point, set it in the icc, */
/* and make the Lut relative to them. */
/* Return NULL on error, check errc+err for reason */
static icxLuBase *xicc_set_luobj(
xicc *p,					/* this */
icmLookupFunc func,			/* Functionality */
icRenderingIntent intent,	/* Intent */
icmLookupOrder order,		/* Search Order */
int flags,					/* white/black point, verbose flags etc. */
int no,						/* Number of points */
int nobw,					/* Number of points to look for white & black patches in */
cow *points,				/* Array of input points in target PCS space */
icxMatrixModel *skm,   		/* Optional skeleton model (used for input profiles) */
double dispLuminance,		/* > 0.0 if display luminance value and is known */
double wpscale,				/* > 0.0 if input white point is to be scaled */
double smooth,				/* RSPL smoothing factor, -ve if raw */
double avgdev,				/* reading Average Deviation as a proportion of the input range */
icxViewCond *vc,			/* Viewing Condition (NULL if not using CAM) */
icxInk *ink,				/* inking details (NULL for default) */
xcal *cal,                  /* Optional cal, will override any existing (not deleted with xicc)*/
int quality					/* Quality metric, 0..3 */
) {
	icmLuBase *plu;
	icxLuBase *xplu = NULL;
	icmLuAlgType alg;

	if (cal != NULL) {
		if (p->cal != NULL && p->nodel_cal == 0)
			p->cal->del(p->cal);
		p->cal = cal;
		p->nodel_cal = 1;	/* We were given it, so don't delete it */
	}

	if (func != icmFwd) {
		p->errc = 1;
		sprintf(p->err,"Can only create Device->PCS profiles from scattered data.");
		xplu = NULL;
		return xplu;
	}

	/* Get icclib lookup object */
	if ((plu = p->pp->get_luobj(p->pp, func, intent, 0, order)) == NULL) {
		p->errc = p->pp->errc;		/* Copy error */
		strcpy(p->err, p->pp->err);
		return NULL;
	}

	/* Figure out what the algorithm is */
	plu->spaces(plu, NULL, NULL, NULL, NULL, &alg, NULL, NULL, NULL, NULL);

	/* Call xiccLu wrapper creation */
	switch (alg) {
		case icmMonoFwdType:
			p->errc = 1;
			sprintf(p->err,"Setting Monochrome Fwd profile from scattered data not supported.");
			plu->del(plu);
			xplu = NULL;		/* Not supported yet */
			break;

    	case icmMatrixFwdType:
			if (smooth < 0.0)
				smooth = -smooth;
			xplu = set_icxLuMatrix(p, plu, flags, no, nobw, points, skm, dispLuminance, wpscale, quality, smooth);
			break;

    	case icmLutType:
			/* ~~~ Should add check that it is a fwd profile ~~~ */
			xplu = set_icxLuLut(p, plu, func, intent, flags, no, nobw, points, skm, dispLuminance, wpscale, smooth, avgdev, vc, ink, quality);
			break;

		default:
			break;
	}

	return xplu;
}

/* ------------------------------------------------------ */
/* Viewing Condition Parameter stuff */

#ifdef NEVER	/* Not currently used */

/* Guess viewing parameters from the technology signature */
static void guess_from_techsig(
icTechnologySignature tsig,
double *Ybp
) {
double Yb = -1.0;

	switch (tsig) {
		/* These are all inputing either a representation of */
		/* a natural scene captured on another medium, or are assuming */
		/* that the medium is the original. A _good_ system would */
		/* let the user indicate which is the case. */
		case icSigReflectiveScanner:
		case icSigFilmScanner:
			Yb = 0.2;
			break;

		/* Direct scene to value devices. */
		case icSigDigitalCamera:
		case icSigVideoCamera:
			Yb = 0.2;
			break;

		/* Emmisive displays. */
		/* We could try tweaking the white point on the assumption */
		/* that the viewer will be adapted to a combination of both */
		/* the CRT white point, and the ambient light. */
		case icSigVideoMonitor:
		case icSigCRTDisplay:
		case icSigPMDisplay:
		case icSigAMDisplay:
			Yb = 0.2;
			break;

		/* Photo CD has its own viewing definitions */
		/* (It represents original scene colors) */
		case icSigPhotoCD:
			Yb = 0.2;
			break;

		/* Projection devices, either direct, or */
		/* via another intermediate medium. */
		case icSigProjectionTelevision:
			Yb = 0.1;		/* Assume darkened room, little background */
			break;
		case icSigFilmWriter:
			Yb = 0.0;		/* Assume a dark room - no background */
			break;

		/* Printed media devices. */
		case icSigInkJetPrinter:
		case icSigThermalWaxPrinter:
		case icSigElectrophotographicPrinter:
		case icSigElectrostaticPrinter:
		case icSigDyeSublimationPrinter:
		case icSigPhotographicPaperPrinter:
		case icSigPhotoImageSetter:
		case icSigGravure:
		case icSigOffsetLithography:
		case icSigSilkscreen:
		case icSigFlexography:
			Yb = 0.2;
			break;

		default:
			Yb = 0.2;
	}

	if (Ybp != NULL)
		*Ybp = Yb;
}

#endif /* NEVER */


/* See if we can read or guess the viewing conditions */
/* for an ICC profile. */
/* Return value 0 if it is well defined */
/* Return value 1 if it is a guess */
/* Return value 2 if it is not possible/appropriate */
int xicc_get_viewcond(
xicc *p,			/* Expanded profile we're working with */
icxViewCond *vc		/* Viewing parameters to return */
) {
	icc *pp = p->pp;	/* Base ICC */

	/* Numbers we're trying to find */
	ViewingCondition Ev = vc_none;
	double Wxyz[3] = {-1.0, -1.0, -1.0};	/* Adapting white color */
	double La = -1.0;		/* Adapting luminance */
	double Ixyz[3] = {-1.0, -1.0, -1.0};	/* Illuminant color */
	double Li = -1.0;						/* Illuminant luminance */
	double Lb = -1.0;		/* Backgrount luminance */
	double Yb = -1.0;		/* Background relative luminance to Lv */
	double Lve = -1.0;		/* Emissive device image luminance */
	double Lvr = -1.0;		/* Reflective device image luminance */
	double Lv = -1.0;		/* device image luminance */
	double Yf = -1.0;		/* Flare relative luminance to Lv */
	double Fxyz[3] = {-1.0, -1.0, -1.0};	/* Flare color */
    icTechnologySignature tsig = icMaxEnumTechnology; /* Technology Signature */
    icProfileClassSignature devc = icMaxEnumClass;
	int trans = -1;		/* Set to 0 if not transparency, 1 if it is */

	/* Collect all the information we can find */

	/* Emmisive devices image white luminance */
	{
		icmXYZArray *luminanceTag;

		if ((luminanceTag = (icmXYZArray *)pp->read_tag(pp, icSigLuminanceTag)) != NULL
         && luminanceTag->ttype == icSigXYZType && luminanceTag->size >= 1) {
			Lve = luminanceTag->data[0].Y;	/* Copy structure */
		} 
	}

	/* Flare: */
	{
		icmMeasurement *ro;

		if ((ro = (icmMeasurement *)pp->read_tag(pp, icSigMeasurementTag)) != NULL
		 && ro->ttype == icSigMeasurementType) {
	
    		Yf = ro->flare;
   			/* ro->illuminant ie D50, D65, D93, A etc. */
		}
	}

	/* Media White Point */
	{
		icmXYZArray *whitePointTag;

		if ((whitePointTag = (icmXYZArray *)pp->read_tag(pp, icSigMediaWhitePointTag)) != NULL
         && whitePointTag->ttype == icSigXYZType && whitePointTag->size >= 1) {
			Wxyz[0] = whitePointTag->data[0].X;
			Wxyz[1] = whitePointTag->data[0].Y;
			Wxyz[2] = whitePointTag->data[0].Z;
		}
	}

	/* ViewingConditions: */
	{
		icmViewingConditions *ro;

		if ((ro = (icmViewingConditions *)pp->read_tag(pp, icSigViewingConditionsTag)) != NULL
		 && ro->ttype == icSigViewingConditionsType) {
	
   			/* ro->illuminant.X */
   			/* ro->illuminant.Z */

    		Li = ro->illuminant.Y;

			/* Reflect illuminant off the media white */
    		Lvr = Li * Wxyz[1];

			/* Illuminant color */
			Ixyz[0] = ro->illuminant.X/ro->illuminant.Y;
			Ixyz[1] = 1.0;
			Ixyz[2] = ro->illuminant.Z/ro->illuminant.Y;
			
			/* Assume ICC surround is CICAM97 background */
    		/* ro->surround.X */
    		/* ro->surround.Z */
    		La = ro->surround.Y;

    		/* ro->stdIlluminant ie D50, D65, D93, A etc. */
		}
	}

	/* Stuff we might need */

	/* Technology: */
	{
		icmSignature *ro;

		/* Try and read the tag from the file */
		if ((ro = (icmSignature *)pp->read_tag(pp, icSigTechnologyTag)) != NULL
		 && ro->ttype != icSigSignatureType) {

		tsig = ro->sig;
		}
	}

    devc = pp->header->deviceClass;	/* Type of profile */
	if (devc == icSigLinkClass
	 || devc == icSigAbstractClass
	 || devc == icSigColorSpaceClass
	 || devc == icSigNamedColorClass)
		return 2;

	/*
    	icSigInputClass
    	icSigDisplayClass
    	icSigOutputClass
	*/

	if ((pp->header->flags & icTransparency) != 0)
		trans = 1;
	else
		trans = 0;


	/* figure Lv if we have the information */
	if (Lve >= 0.0)
		Lv = Lve;		/* Emmisive image white luminance */
	else
		Lv = Lvr;		/* Reflectance image white luminance */

	/* Fudge the technology signature */
	if (tsig == icMaxEnumTechnology) {
		if (devc == icSigDisplayClass)
			tsig = icSigCRTDisplay;
	}

#ifndef NEVER
	printf("Enumeration = %d\n", Ev);
	printf("Viewing Conditions:\n");
	printf("White adaptation color %f %f %f\n",Wxyz[0], Wxyz[1], Wxyz[2]);
	printf("Adapting Luminance La = %f\n",La);
	printf("Illuminant color %f %f %f\n",Ixyz[0], Ixyz[1], Ixyz[2]);
	printf("Illuminant Luminance Li = %f\n",Li);
	printf("Background Luminance Lb = %f\n",Lb);
	printf("Relative Background Yb = %f\n",Yb);
	printf("Emissive Image White Lve = %f\n",Lve);
	printf("Reflective Image White Lvr = %f\n",Lvr);
	printf("Device Image White Lv = %f\n",Lv);
	printf("Relative Flare Yf = %f\n",Yf);
	printf("Flare color %f %f %f\n",Fxyz[0], Fxyz[1], Fxyz[2]);
	printf("Technology = %s\n",tag2str(tsig));
	printf("deviceClass = %s\n",tag2str(devc));
	printf("Transparency = %d\n",trans);
#endif
	
	/* See if the viewing conditions are completely defined as ICC can do it */
	if (Wxyz[0] >= 0.0 && Wxyz[1] >= 0.0 && Wxyz[2] >= 0.0
	 && La >= 0.0
	 && Yb >= 0.0
	 && Lv >= 0.0
	 && Yf >= 0.0
	 && Fxyz[0] >= 0.0 && Fxyz[1] >= 0.0 && Fxyz[2] >= 0.0) {

		vc->Ev = vc_none;
		vc->Wxyz[0] = Wxyz[0];
		vc->Wxyz[1] = Wxyz[1];
		vc->Wxyz[2] = Wxyz[2];
		vc->La = La;
		vc->Yb = Yb;
		vc->Lv = Lv;
		vc->Yf = Yf;
		vc->Fxyz[0] = Fxyz[0];	
		vc->Fxyz[1] = Fxyz[1];	
		vc->Fxyz[2] = Fxyz[2];	
		return 0;
	}

	/* Hmm. We didn't get all the info an ICC can contain. */
	/* We will try to guess some reasonable defaults */

	/* Have we at least got an adaptation white point ? */
	if (Wxyz[0] < 0.0 || Wxyz[1] < 0.0 || Wxyz[2] < 0.0)
		return 2;						/* No */

	/* Have we got the technology ? */
	if (tsig == icMaxEnumTechnology)
		return 2;						/* Hopeless */

	/* Guess from the technology */
	switch (tsig) {

		/* This is inputing either a representation of */
		/* a natural scene captured on another a print medium, or */
		/* are is assuming that the medium is the original. */
		/* We will assume that the print is the original. */
		case icSigReflectiveScanner:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 34.0;	/* Use a practical print evaluation number */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.2;	/* Assume grey world */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_average;	/* Assume average viewing conditions */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.01;	/* Assume 1% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}

		/* This is inputing either a representation of */
		/* a natural scene captured on another a photo medium, or */
		/* are is assuming that the medium is the original. */
		/* We will assume a compromise media original, natural scene */
		case icSigFilmScanner:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 50.0;	/* Use bright indoors, dull outdoors */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.2;	/* Assume grey world */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_average;	/* Assume average viewing conditions */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.005;	/* Assume 0.5% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}

		/* Direct scene to value devices. */
		case icSigDigitalCamera:
		case icSigVideoCamera:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 110.0;	/* Use very bright indoors, usual outdoors */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.2;	/* Assume grey world */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_average;	/* Assume average viewing conditions */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.0;	/* Assume 0% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}

		/* Emmisive displays. */
		/* Assume a video monitor is in a darker environment than a CRT */
		case icSigVideoMonitor:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 4.0;	/* Darkened work environment */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.2;	/* Assume grey world */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_dim;	/* Assume dim viewing conditions */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.01;	/* Assume 1% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}


		/* Assume a typical work environment */
		case icSigCRTDisplay:
		case icSigPMDisplay:
		case icSigAMDisplay:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 33.0;	/* Typical work environment */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.2;	/* Assume grey world */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_average;	/* Assume average viewing conditions */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.02;	/* Assume 2% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}

		/* Photo CD has its own viewing definitions */
		/* (It represents original scene colors) */
		case icSigPhotoCD:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 320.0;	/* Bright outdoors */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.2;	/* Assume grey world */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_average;	/* Assume average viewing conditions */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.00;	/* Assume 0% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}

		/* Projection devices, either direct, or */
		/* via another intermediate medium. */
		/* Assume darkened room, little background */
		case icSigProjectionTelevision:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 7.0;	/* Dark environment */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.1;	/* Assume little background */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_dim;	/* Dim environment */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.01;	/* Assume 1% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}
		/* Assume very darkened room, no background */
		case icSigFilmWriter:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 7.0;	/* Dark environment */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.0;	/* Assume no background */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_dark;	/* Dark environment */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.01;	/* Assume 1% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}

		/* Printed media devices. */
		/* Assume a normal print viewing environment */
		case icSigInkJetPrinter:
		case icSigThermalWaxPrinter:
		case icSigElectrophotographicPrinter:
		case icSigElectrostaticPrinter:
		case icSigDyeSublimationPrinter:
		case icSigPhotographicPaperPrinter:
		case icSigPhotoImageSetter:
		case icSigGravure:
		case icSigOffsetLithography:
		case icSigSilkscreen:
		case icSigFlexography:
			{
				if (La < 0.0)	/* No adapting luminance */
					La = 40.0;	/* Use a practical print evaluation number */
				if (Yb < 0.0)	/* No background relative luminance */
					Yb = 0.2;	/* Assume grey world */
				if (Lv < 0.0)	/* No device image luminance */
					Ev = vc_average;	/* Assume average viewing conditions */
				if (Yf < 0.0)	/* No flare figure */
					Yf = 0.01;	/* Assume 1% flare */
				if (Fxyz[0] < 0.0 || Fxyz[1] < 0.0 || Fxyz[2] < 0.0)	/* No flare color */
					Fxyz[0] = Wxyz[0], Fxyz[1] = Wxyz[1], Fxyz[2] = Wxyz[2];
				break;
			}

		default:
			{
			return 2;
			}
	}

	return 1;
}

/* Write our viewing conditions to the underlying ICC profile, */
/* using a private tag. */
void xicc_set_viewcond(
xicc *p,			/* Expanded profile we're working with */
icxViewCond *vc		/* Viewing parameters to return */
) {
	//icc *pp = p->pp;	/* Base ICC */

	// ~~1	Not implemented yet
}



/* Return an enumerated viewing condition */
/* Return enumeration if OK, -999 if there is no such enumeration. */
/* xicc may be NULL if just the description is wanted, */
/* or an explicit white point is provided. */
int xicc_enum_viewcond(
xicc *p,			/* Expanded profile to get white point (May be NULL if desc NZ) */
icxViewCond *vc,	/* Viewing parameters to return, May be NULL if desc is nz */
int no,				/* Enumeration to return, -1 for default, -2 for none */
char *as,			/* String alias to number, NULL if none */
int desc,			/* NZ - Just return a description of this enumeration in vc */
double *wp			/* Provide white point if xicc is NULL */
) {

	if (desc == 0) {	/* We're setting the viewing condition */
		icc *pp;		/* Base ICC */
		icmXYZArray *whitePointTag;

		if (vc == NULL)
			return -999;

		if (p == NULL) {
			if (wp == NULL)
				return -999;
			vc->Wxyz[0] = wp[0];
			vc->Wxyz[1] = wp[1];
			vc->Wxyz[2] = wp[2];
		} else {
	
			pp = p->pp;
			if ((whitePointTag = (icmXYZArray *)pp->read_tag(pp, icSigMediaWhitePointTag)) != NULL
	   	      && whitePointTag->ttype == icSigXYZType && whitePointTag->size >= 1) {
				vc->Wxyz[0] = whitePointTag->data[0].X;
				vc->Wxyz[1] = whitePointTag->data[0].Y;
				vc->Wxyz[2] = whitePointTag->data[0].Z;
			} else {
				if (wp == NULL) { 
					sprintf(p->err,"Enum VC: Failed to read Media White point");
					p->errc = 2;
					return -999;
				}
				vc->Wxyz[0] = wp[0];
				vc->Wxyz[1] = wp[1];
				vc->Wxyz[2] = wp[2];
			}
		}

		/* Set a default flare color */
		vc->Fxyz[0] = vc->Wxyz[0];
		vc->Fxyz[1] = vc->Wxyz[1];
		vc->Fxyz[2] = vc->Wxyz[2];
	}

	/*

	Typical adapting field luminances and white luminance in reflective setup:

	E = illuminance in Lux
	Lv = White luminance assuming 100% reflectance
	La = Adapting field luminance in cd/m^2, assuming 20% reflectance from surround
    
	    E     La     Lv   Condition 
	    11     0.7    4   Twilight
	    32     2     10   Subdued indoor lighting
	    64     4     20   Less than typical office light; sometimes recommended for
	                      display-only workplaces (sRGB)
	   350    22    111   Typical Office (sRGB annex D)
	   500    32    160   Practical print evaluationa (ISO-3664 P2)
	  1000    64    318   Good Print evaluation (CIE 116-1995)
	  1000    64    318   Television Studio lighting
	  1000    64    318   Overcast Outdoors
	  2000   127    637   Critical print evaluation (ISO-3664 P1)
	 10000   637   3183   Typical outdoors, full daylight 
	 50000  3185  15915   Bright summers day 

	*/

	if (no == -1
	 || (as != NULL && stricmp(as,"d") == 0)) {

		no = -1;
		if (vc != NULL) {
			vc->desc = "  d - Default Viewing Condition";
			vc->Ev = vc_average;	/* Average viewing conditions */
			vc->La = 50.0;			/* Practical to Good lighting */
			vc->Lv = 250.0;			/* Average viewing conditions ratio */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.01;			/* 1% flare */
		}
	}
	else if (no == 0
	 || (as != NULL && stricmp(as,"pp") == 0)) {

		no = 0;
		if (vc != NULL) {
			vc->desc = " pp - Practical Reflection Print (ISO-3664 P2)";
			vc->Ev = vc_average;	/* Average viewing conditions */
			vc->La = 32.0;			/* Use a practical print evaluation number */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.01;			/* 1% flare */
		}
	}
	else if (no == 1
	 || (as != NULL && stricmp(as,"pe") == 0)) {

		no = 1;
		if (vc != NULL) {
			vc->desc = " pe - Print evaluation environment (CIE 116-1995)";
			vc->Ev = vc_average;	/* Average viewing conditions */
			vc->La = 64.0;			/* Good */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.01;			/* 1% flare */
		}
	}
	else if (no == 2
	 || (as != NULL && stricmp(as,"pc") == 0)) {

		no = 2;
		if (vc != NULL) {
			vc->desc = " pc - Critical print evaluation environment (ISO-3664 P1)";
			vc->Ev = vc_average;	/* Average viewing conditions */
			vc->La = 127.0;			/* Critical */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.01;			/* 1% flare */
		}
	}
	else if (no == 3
	 || (as != NULL && stricmp(as,"mt") == 0)) {

		no = 3;
		if (vc != NULL) {
			vc->desc = " mt - Monitor in typical work environment";
			vc->Ev = vc_average;	/* Average viewing conditions */
			vc->La = 22.0;			/* Typical work environment */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.02;			/* 2% flare */
		}
	}
	else if (no == 4
	 || (as != NULL && stricmp(as,"mb") == 0)) {

		no = 4;
		if (vc != NULL) {
			vc->desc = " mb - Bright monitor in bright work environment";
			vc->Ev = vc_average;	/* Average viewing conditions */
			vc->La = 42.0;			/* Bright work environment */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.02;			/* 2% flare */
		}
	}
	else if (no == 5
	 || (as != NULL && stricmp(as,"md") == 0)) {

		no = 5;
		if (vc != NULL) {
			vc->desc = " md - Monitor in darkened work environment";
			vc->Ev = vc_dim;		/* Dim viewing conditions */
			vc->La = 4.0;			/* Darkened work environment */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.01;			/* 1% flare */
		}
	}
	else if (no == 6
	 || (as != NULL && stricmp(as,"jm") == 0)) {

		no = 6;
		if (vc != NULL) {
			vc->desc = " jm - Projector in dim environment";
			vc->Ev = vc_dim;		/* Dim viewing conditions */
			vc->La = 10.0;			/* Adaptation is from display */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.01;			/* 1% flare */
		}
	}
	else if (no == 7
	 || (as != NULL && stricmp(as,"jd") == 0)) {

		no = 7;
		if (vc != NULL) {
			vc->desc = " jd - Projector in dark environment";
			vc->Ev = vc_dark;		/* Dark viewing conditions */
			vc->La = 10.0;			/* Adaptation is from display */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.01;			/* 1% flare ? */
		}
	}
	else if (no == 8
	 || (as != NULL && stricmp(as,"pcd") == 0)) {

		no = 8;
		if (vc != NULL) {
			vc->desc = "pcd - Photo CD - original scene outdoors";
			vc->Ev = vc_average;	/* Average viewing conditions */
			vc->La = 320.0;			/* Typical outdoors, 1600 cd/m^2 */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.00;			/* 0% flare */
		}
	}
	else if (no == 9
	 || (as != NULL && stricmp(as,"ob") == 0)) {

		no = 9;
		if (vc != NULL) {
			vc->desc = " ob - Original scene - Bright Outdoors";
			vc->Ev = vc_average;	/* Average viewing conditions */
			vc->La = 2000.0;		/* Bright Outdoors */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.00;			/* 0% flare */
		}
	}
	else if (no == 10
	 || (as != NULL && stricmp(as,"cx") == 0)) {

		no = 10;
		if (vc != NULL) {
			vc->desc = " cx - Cut Sheet Transparencies on a viewing box";
			vc->Ev = vc_cut_sheet;	/* Cut sheet viewing conditions */
			vc->La = 53.0;			/* Dim, adapted to slide ? */
			vc->Yb = 0.2;			/* Grey world */
			vc->Yf = 0.01;			/* 1% flare ? */
		}
	}
	else {
		if (p != NULL) {
			sprintf(p->err,"Enum VC: Unrecognised enumeration %d",no);
			p->errc = 1;
		}
		return -999;
	}

	return no;
}

/* Debug: dump a Viewing Condition to standard out */
void xicc_dump_viewcond(
icxViewCond *vc
) {
	printf("Viewing Condition:\n");
	if (vc->Ev == vc_dark)
		printf("  Surround to Image: Dark\n");
	else if (vc->Ev == vc_dim)
		printf("  Surround to Image: Dim\n");
	else if (vc->Ev == vc_average)
		printf("  Surround to Image: Average\n");
	else if (vc->Ev == vc_cut_sheet)
		printf("  Transparency on Light box\n");
	printf("  Adapted white = %f %f %f\n",vc->Wxyz[0], vc->Wxyz[1], vc->Wxyz[2]);
	printf("  Adapted luminance = %f cd/m^2\n",vc->La);
	printf("  Background to image ratio = %f\n",vc->Yb);
	if (vc->Ev == vc_none)
		printf("  Image luminance = %f cd/m^2\n",vc->Lv);
	printf("  Flare to image ratio = %f\n",vc->Yf);
	printf("  Flare color = %f %f %f\n",vc->Fxyz[0], vc->Fxyz[1], vc->Fxyz[2]);
}


/* Debug: dump an Inking setup to standard out */
void xicc_dump_inking(icxInk *ik) {
	printf("Inking settings:\n");
	if (ik->tlimit < 0.0)
		printf("No total limit\n");
	else
		printf("Total limit = %f%%\n",ik->tlimit * 100.0);

	if (ik->klimit < 0.0)
		printf("No black limit\n");
	else
		printf("Black limit = %f%%\n",ik->klimit * 100.0);

	if (ik->KonlyLmin)
		printf("K only black as locus Lmin\n");
	else
		printf("Normal black as locus Lmin\n");

	if (ik->k_rule == icxKvalue) {
		printf("Inking rule is a fixed K target\n");
	} if (ik->k_rule == icxKlocus) {
		printf("Inking rule is a fixed locus target\n");
	} if (ik->k_rule == icxKluma5 || ik->k_rule == icxKluma5k) {
		if (ik->k_rule == icxKluma5)
			printf("Inking rule is a 5 parameter locus function of L\n");
		else
			printf("Inking rule is a 5 parameter K function of L\n");
		printf("Ksmth = %f\n",ik->c.Ksmth);
		printf("Kskew = %f\n",ik->c.Kskew);
		printf("Kstle = %f\n",ik->c.Kstle);
		printf("Kstpo = %f\n",ik->c.Kstpo);
		printf("Kenpo = %f\n",ik->c.Kenpo);
		printf("Kenle = %f\n",ik->c.Kenle);
		printf("Kshap = %f\n",ik->c.Kshap);
	} if (ik->k_rule == icxKl5l || ik->k_rule == icxKl5lk) {
		if (ik->k_rule == icxKl5l)
			printf("Inking rule is a 2x5 parameter locus function of L and K aux\n");
		else
			printf("Inking rule is a 2x5 parameter K function of L and K aux\n");
		printf("Min Ksmth = %f\n",ik->c.Ksmth);
		printf("Min Kskew = %f\n",ik->c.Kskew);
		printf("Min Kstle = %f\n",ik->c.Kstle);
		printf("Min Kstpo = %f\n",ik->c.Kstpo);
		printf("Min Kenpo = %f\n",ik->c.Kenpo);
		printf("Min Kenle = %f\n",ik->c.Kenle);
		printf("Min Kshap = %f\n",ik->c.Kshap);
		printf("Max Ksmth = %f\n",ik->x.Ksmth);
		printf("Max Kskew = %f\n",ik->x.Kskew);
		printf("Max Kstle = %f\n",ik->x.Kstle);
		printf("Max Kstpo = %f\n",ik->x.Kstpo);
		printf("Max Kenpo = %f\n",ik->x.Kenpo);
		printf("Max Kenle = %f\n",ik->x.Kenle);
		printf("Max Kshap = %f\n",ik->x.Kshap);
	} 
}

/* ------------------------------------------------------ */
/* Gamut Mapping Intent stuff */

/* Return an enumerated gamut mapping intent */
/* Return enumeration if OK, icxIllegalGMIntent if there is no such enumeration. */
int xicc_enum_gmapintent(
icxGMappingIntent *gmi,	/* Gamut Mapping parameters to return */
int no,					/* Enumeration selected, icxNoGMIntent for none */
char *as				/* Alias string selector, NULL for none */
) {
#ifdef USE_CAM
	int colccas = 0x2;	/* Use cas clipping for colorimetric style intents */
	int perccas = 0x1;	/* Use cas for perceptual style intents */
#else
	int colccas = 0x0;	/* Use Lab for colorimetric style intents */
	int perccas = 0x0;	/* Use Lab for perceptual style intents */
	fprintf(stderr,"!!!!!! Warning, USE_CAM is off in xicc.c !!!!!!\n");
#endif

	/* Assert default if no guidance given */
	if (no == icxNoGMIntent && as == NULL)
		no = icxDefaultGMIntent;

	if (no == 0
	 || no == icxAbsoluteGMIntent
	 || (as != NULL && stricmp(as,"a") == 0)) {
		/* Map Absolute appearance space Jab to Jab and clip out of gamut */
		no = 0;
		gmi->as = "a";
		gmi->desc = " a - Absolute Colorimetric (in Jab) [ICC Absolute Colorimetric]";
		gmi->icci = icAbsoluteColorimetric;
		gmi->usecas  = colccas;		/* Use absolute appearance space */
		gmi->usemap  = 0;			/* Don't use gamut mapping */
		gmi->greymf  = 0.0;
		gmi->glumwcpf = 0.0;
		gmi->glumwexf = 0.0;
		gmi->glumbcpf = 0.0;
		gmi->glumbexf = 0.0;
		gmi->glumknf = 0.0;
		gmi->gamcpf  = 0.0;
		gmi->gamexf  = 0.0;
		gmi->gamcknf  = 0.0;
		gmi->gamxknf  = 0.0;
		gmi->gampwf  = 0.0;
		gmi->gamswf  = 0.0;
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 1
	 || (as != NULL && stricmp(as,"aw") == 0)) {

		/* I'm not sure how often this intent is useful. It's less likely than */
		/* I though that a printer white point won't fit within the gamut */
		/* of a display profile, since the display white always has Y = 1.0, */
		/* and no paper has better than about 95% reflectance. */
		/* Perhaps it may be more useful for targeting printer profiles ? */

		/* Map Absolute Jab to Jab and scale source to avoid clipping the white point */
		no = 1;
		gmi->as = "aw";
		gmi->desc = "aw - Absolute Colorimetric (in Jab) with scaling to fit white point";
		gmi->icci = icAbsoluteColorimetric;
		gmi->usecas = 0x100 | colccas;	/* Absolute Appearance space with scaling */
									/* to avoid clipping the source white point */
		gmi->usemap  = 0;			/* Don't use gamut mapping */
		gmi->greymf  = 0.0;
		gmi->glumwcpf = 0.0;
		gmi->glumwexf = 0.0;
		gmi->glumbcpf = 0.0;
		gmi->glumbexf = 0.0;
		gmi->glumknf = 0.0;
		gmi->gamcpf  = 0.0;
		gmi->gamexf  = 0.0;
		gmi->gamcknf  = 0.0;
		gmi->gamxknf  = 0.0;
		gmi->gampwf  = 0.0;
		gmi->gamswf  = 0.0;
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 2
	 || (as != NULL && stricmp(as,"aa") == 0)) {

		/* Map appearance space Jab to Jab and clip out of gamut */
		no = 2;
		gmi->as = "aa";
		gmi->desc = "aa - Absolute Appearance";
		gmi->icci = icRelativeColorimetric;
		gmi->usecas  = perccas;		/* Appearance space */
		gmi->usemap  = 0;			/* Don't use gamut mapping */
		gmi->greymf  = 0.0;
		gmi->glumwcpf = 0.0;
		gmi->glumwexf = 0.0;
		gmi->glumbcpf = 0.0;
		gmi->glumbexf = 0.0;
		gmi->glumknf = 0.0;
		gmi->gamcpf  = 0.0;
		gmi->gamexf  = 0.0;
		gmi->gamcknf  = 0.0;
		gmi->gamxknf  = 0.0;
		gmi->gampwf  = 0.0;
		gmi->gamswf  = 0.0;
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 3
	 || no == icxRelativeGMIntent
	 || (as != NULL && stricmp(as,"r") == 0)) {

		/* Align neutral axes and linearly map white point, then */
		/* map appearance space Jab to Jab and clip out of gamut */
		no = 3;
		gmi->as = "r";
		gmi->desc = " r - White Point Matched Appearance [ICC Relative Colorimetric]";
		gmi->icci = icRelativeColorimetric;
		gmi->usecas  = perccas;		/* Appearance space */
		gmi->usemap  = 1;			/* Use gamut mapping */
		gmi->greymf  = 1.0;			/* Fully align grey axis */
		gmi->glumwcpf = 1.0;		/* Fully compress grey axis at white end */
		gmi->glumwexf = 1.0;		/* Fully expand grey axis at white end */
		gmi->glumbcpf = 0.0;		/* No compression at black end */
		gmi->glumbexf = 0.0;		/* No expansion at black end */
		gmi->glumknf = 0.0;
		gmi->gamcpf  = 0.0;
		gmi->gamexf  = 0.0;
		gmi->gamcknf  = 0.0;
		gmi->gamxknf  = 0.0;
		gmi->gampwf  = 0.0;
		gmi->gamswf  = 0.0;
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 4
	 || (as != NULL && stricmp(as,"la") == 0)) {

		/* Align neutral axes and linearly map white and black points, then */
		/* map appearance space Jab to Jab and clip out of gamut */
		no = 4;
		gmi->as = "la";
		gmi->desc = "la - Luminance axis matched Appearance";
		gmi->icci = icRelativeColorimetric;
		gmi->usecas  = perccas;		/* Appearance space */
		gmi->usemap  = 1;			/* Use gamut mapping */
		gmi->greymf  = 1.0;			/* Fully align grey axis */
		gmi->glumwcpf = 1.0;		/* Fully compress grey axis at white end */
		gmi->glumwexf = 1.0;		/* Fully expand grey axis at white end */
		gmi->glumbcpf = 1.0;		/* Fully compress grey axis at black end */
		gmi->glumbexf = 1.0;		/* Fully expand grey axis at black end */
		gmi->glumknf = 0.0;			/* No knee on grey mapping */
		gmi->gamcpf  = 0.0;			/* No gamut compression */
		gmi->gamexf  = 0.0;			/* No gamut expansion */
		gmi->gamcknf  = 0.0;		/* No knee in gamut compress */
		gmi->gamxknf  = 0.0;		/* No knee in gamut expand */
		gmi->gampwf  = 0.0;			/* No Perceptual surface weighting factor */
		gmi->gamswf  = 0.0;			/* No Saturation surface weighting factor */
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 5
	 || no == icxDefaultGMIntent
	 || no == icxPerceptualGMIntent
	 || (as != NULL && stricmp(as,"p") == 0)) {

		/* Align neutral axes and perceptually map white and black points, */ 
		/* perceptually compress out of gamut and map appearance space Jab to Jab. */  
		no = 5;
		gmi->as = "p";
		gmi->desc = " p - Perceptual (Preferred) (Default) [ICC Perceptual]";
		gmi->icci = icPerceptual;
		gmi->usecas  = perccas;		/* Appearance space */
		gmi->usemap  = 1;			/* Use gamut mapping */
		gmi->greymf  = 1.0;			/* Fully align grey axis */
		gmi->glumwcpf = 1.0;		/* Fully compress grey axis at white end */
		gmi->glumwexf = 1.0;		/* Fully expand grey axis at white end */
		gmi->glumbcpf = 1.0;		/* Fully compress grey axis at black end */
		gmi->glumbexf = 1.0;		/* Fully expand grey axis at black end */
		gmi->glumknf = 1.0;			/* Sigma knee in grey compress/expand */
		gmi->gamcpf  = 1.0;			/* Full gamut compression */
		gmi->gamexf  = 0.0;			/* No gamut expansion */
		gmi->gamcknf  = 0.8;		/* High Sigma knee in gamut compress */
		gmi->gamxknf  = 0.0;		/* No knee in gamut expand */
		gmi->gampwf  = 1.0;			/* Full Perceptual surface weighting factor */
		gmi->gamswf  = 0.0;			/* No Saturation surface weighting factor */
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 6
	 || (as != NULL && stricmp(as,"pa") == 0)) {

		/* Don't align neutral axes, but perceptually compress out of gamut */
		/* and map appearance space Jab to Jab. */  
		no = 5;
		gmi->as = "pa";
		gmi->desc = "pa - Perceptual Apperance ";
		gmi->icci = icPerceptual;
		gmi->usecas  = perccas;		/* Appearance space */
		gmi->usemap  = 1;			/* Use gamut mapping */
		gmi->greymf  = 0.0;			/* Don't align grey axis */
		gmi->glumwcpf = 1.0;		/* Fully compress grey axis at white end */
		gmi->glumwexf = 1.0;		/* Fully expand grey axis at white end */
		gmi->glumbcpf = 1.0;		/* Fully compress grey axis at black end */
		gmi->glumbexf = 1.0;		/* Fully expand grey axis at black end */
		gmi->glumknf = 1.0;			/* Sigma knee in grey compress/expand */
		gmi->gamcpf  = 1.0;			/* Full gamut compression */
		gmi->gamexf  = 0.0;			/* No gamut expansion */
		gmi->gamcknf  = 0.8;		/* High Sigma knee in gamut compress */
		gmi->gamxknf  = 0.0;		/* No knee in gamut expand */
		gmi->gampwf  = 1.0;			/* Full Perceptual surface weighting factor */
		gmi->gamswf  = 0.0;			/* No Saturation surface weighting factor */
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 7
	 || (as != NULL && stricmp(as,"ms") == 0)) {

		/* Align neutral axes and perceptually map white and black points, */ 
		/* perceptually compress and expand to match gamuts and map Jab to Jab. */  
		no = 6;
		gmi->as = "ms";
		gmi->desc = "ms - Saturation";
		gmi->icci = icSaturation;
		gmi->usecas  = perccas;		/* Appearance space */
		gmi->usemap  = 1;			/* Use gamut mapping */
		gmi->greymf  = 1.0;			/* Fully align grey axis */
		gmi->glumwcpf = 1.0;		/* Fully compress grey axis at white end */
		gmi->glumwexf = 1.0;		/* Fully expand grey axis at white end */
		gmi->glumbcpf = 1.0;		/* Fully compress grey axis at black end */
		gmi->glumbexf = 1.0;		/* Fully expand grey axis at black end */
		gmi->glumknf = 1.0;			/* Sigma knee in grey compress/expand */
		gmi->gamcpf  = 1.0;			/* Full gamut compression */
		gmi->gamexf  = 1.0;			/* Full gamut expansion  */
		gmi->gamcknf = 1.0;			/* High Sigma knee in gamut compress/expand */
		gmi->gamxknf = 0.4;			/* Moderate Sigma knee in gamut compress/expand */
		gmi->gampwf  = 0.2;			/* Slight perceptual surface weighting factor */
		gmi->gamswf  = 0.8;			/* Most saturation surface weighting factor */
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 8
	 || no == icxSaturationGMIntent
	 || (as != NULL && stricmp(as,"s") == 0)) {

		/* Same as "ms" but enhance saturation */
		no = 7;
		gmi->as = "s";
		gmi->desc = " s - Enhanced Saturation [ICC Saturation]";
		gmi->icci = icSaturation;
		gmi->usecas  = perccas;		/* Appearance space */
		gmi->usemap  = 1;			/* Use gamut mapping */
		gmi->greymf  = 1.0;			/* Fully align grey axis */
		gmi->glumwcpf = 1.0;		/* Fully compress grey axis at white end */
		gmi->glumwexf = 1.0;		/* Fully expand grey axis at white end */
		gmi->glumbcpf = 1.0;		/* Fully compress grey axis at black end */
		gmi->glumbexf = 1.0;		/* Fully expand grey axis at black end */
		gmi->glumknf = 1.0;			/* Sigma knee in grey compress/expand */
		gmi->gamcpf  = 1.0;			/* Full gamut compression */
		gmi->gamexf  = 1.0;			/* Full gamut expansion */
		gmi->gamcknf = 1.0;			/* High sigma knee in gamut compress */
		gmi->gamxknf = 0.5;			/* Moderate sigma knee in gamut expand */
		gmi->gampwf  = 0.0;			/* No Perceptual surface weighting factor */
		gmi->gamswf  = 1.0;			/* Full Saturation surface weighting factor */
		gmi->satenh  = 0.9;			/* Medium saturation enhancement */
	}
	else if (no == 9
	 || (as != NULL && stricmp(as,"al") == 0)) {

		/* Map absolute L*a*b* to L*a*b* and clip out of gamut */
		no = 8;
		gmi->as = "al";
		gmi->desc = "al - Absolute Colorimetric (Lab)";
		gmi->icci = icAbsoluteColorimetric;
		gmi->usecas  = 0x0;			/* Don't use appearance space, use L*a*b* */
		gmi->usemap  = 0;			/* Don't use gamut mapping */
		gmi->greymf  = 0.0;
		gmi->glumwcpf = 0.0;
		gmi->glumwexf = 0.0;
		gmi->glumbcpf = 0.0;
		gmi->glumbexf = 0.0;
		gmi->glumknf = 0.0;
		gmi->gamcpf  = 0.0;
		gmi->gamexf  = 0.0;
		gmi->gamcknf = 0.0;
		gmi->gamxknf = 0.0;
		gmi->gampwf  = 0.0;
		gmi->gamswf  = 0.0;
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else if (no == 10
	 || (as != NULL && stricmp(as,"rl") == 0)) {

		/* Align neutral axes and linearly map white point, then */
		/* map L*a*b* to L*a*b* and clip out of gamut */
		no = 3;
		gmi->as = "rl";
		gmi->desc = "rl - White Point Matched Appearance (Lab)";
		gmi->icci = icRelativeColorimetric;
		gmi->usecas  = 0x0;			/* Don't use appearance space, use L*a*b* */
		gmi->usemap  = 1;			/* Use gamut mapping */
		gmi->greymf  = 1.0;			/* And linearly map white point */
		gmi->glumwcpf = 1.0;
		gmi->glumwexf = 1.0;
		gmi->glumbcpf = 0.0;
		gmi->glumbexf = 0.0;
		gmi->glumknf = 0.0;
		gmi->gamcpf  = 0.0;
		gmi->gamexf  = 0.0;
		gmi->gamcknf = 0.0;
		gmi->gamxknf = 0.0;
		gmi->gampwf  = 0.0;
		gmi->gamswf  = 0.0;
		gmi->satenh  = 0.0;			/* No saturation enhancement */
	}
	else {	/* icxIllegalGMIntent */
		return icxIllegalGMIntent;
	}

	return no;
}


/* Debug: dump a Gamut Mapping specification */
void xicc_dump_gmi(
icxGMappingIntent *gmi	/* Gamut Mapping parameters to return */
) {
	printf(" Gamut Mapping Specification:\n");
	if (gmi->desc != NULL)
		printf("  Description = '%s'\n",gmi->desc);
	printf("  Closest ICC intent = '%s'\n",icm2str(icmRenderingIntent,gmi->icci));

	if ((gmi->usecas & 0xff) == 0)
		printf("  Not using Color Apperance Space\n");
	else if ((gmi->usecas & 0xff) == 1)
		printf("  Using Color Apperance Space\n");
	else if ((gmi->usecas & 0xff) == 2)
		printf("  Using Absolute Color Apperance Space\n");

	if ((gmi->usecas & 0x100) != 0)
		printf("  Scaling source to avoid white point clipping\n");

	if (gmi->usemap == 0)
		printf("  Not using Mapping\n");
	else {
		printf("  Using Mapping with parameters:\n");
		printf("  Grey axis alignment   factor %f\n", gmi->greymf);
		printf("  Grey axis white compression factor %f\n", gmi->glumwcpf);
		printf("  Grey axis white expansion   factor %f\n", gmi->glumwexf);
		printf("  Grey axis black compression factor %f\n", gmi->glumbcpf);
		printf("  Grey axis black expansion   factor %f\n", gmi->glumbexf);
		printf("  Grey axis knee        factor %f\n", gmi->glumknf);
		printf("  Gamut compression factor %f\n", gmi->gamcpf);
		printf("  Gamut expansion   factor %f\n", gmi->gamexf);
		printf("  Gamut compression knee factor %f\n", gmi->gamcknf);
		printf("  Gamut expansion   knee factor %f\n", gmi->gamxknf);
		printf("  Gamut Perceptual mapping weighting factor %f\n", gmi->gampwf);
		printf("  Gamut Saturation mapping weighting factor %f\n", gmi->gamswf);
		printf("  Saturation enhancement factor %f\n", gmi->satenh);
	}
}

/* ------------------------------------------------------ */
/* Turn xicc xcal into limit calibration callback */

/* Given an icc profile, try and create an xcal */
/* Return NULL on error or no cal */
xcal *xiccReadCalTag(icc *p) {
	xcal *cal = NULL;
	icTagSignature sig = icmMakeTag('t','a','r','g');
	icmText *ro;
	int oi, tab;

//printf("~1 about to look for CAL in profile\n");
	if ((ro = (icmText *)p->read_tag(p, sig)) != NULL) { 
		cgatsFile *cgf;
		cgats *icg;

		if (ro->ttype != icSigTextType)
			return NULL;
	
//printf("~1 found 'targ' tag\n");
		if ((icg = new_cgats()) == NULL) {
			return NULL;
		}
		if ((cgf = new_cgatsFileMem(ro->data, ro->size)) != NULL) {
			icg->add_other(icg, "CTI3");
			oi = icg->add_other(icg, "CAL");

//printf("~1 created cgats object from 'targ' tag\n");
			if (icg->read(icg, cgf) == 0) {

				for (tab = 0; tab < icg->ntables; tab++) {
					if (icg->t[tab].tt == tt_other && icg->t[tab].oi == oi) {
						break;
					}
				}
				if (tab < icg->ntables) {
//printf("~1 found CAL table\n");
		
					if ((cal = new_xcal()) == NULL) {
						icg->del(icg);
						cgf->del(cgf);
						return NULL;
					}
					if (cal->read_cgats(cal, icg, tab, "'targ' tag") != 0)  {
#ifdef DEBUG
						printf("read_cgats on cal tag failed\n");
#endif
						cal->del(cal);
						cal = NULL;
					}
//else printf("~1 read CAL and creaded xcal object OK\n");
				}
			}
			cgf->del(cgf);
		}
		icg->del(icg);
	}
	return cal;
}

/* A callback that uses an xcal, that can be used with icc get_tac */
void xiccCalCallback(void *cntx, double *out, double *in) {
	xcal *cal = (xcal *)cntx;

	cal->interp(cal, out, in);
}

/* ---------------------------------------------- */

/* Utility function - given an open icc profile, */
/* guess which channel is the black. */
/* Return -1 if there is no black channel or it can't be guessed */
int icxGuessBlackChan(icc *p) {
	int kch = -1;

	switch (p->header->colorSpace) {
   	 	case icSigCmykData:
			kch = 3;
			break;

		/* Use a heuristic to detect the black channel. */
		/* This duplicates code in icxLu_comp_bk_point() :-( */ 
		/* Colorant guessing should go in icclib ? */
		case icSig2colorData:
		case icSig3colorData:
		case icSig4colorData:
		case icSig5colorData:
		case icSig6colorData:
		case icSig7colorData:
		case icSig8colorData:
		case icSig9colorData:
		case icSig10colorData:
		case icSig11colorData:
		case icSig12colorData:
		case icSig13colorData:
		case icSig14colorData:
		case icSig15colorData:
		case icSigMch5Data:
		case icSigMch6Data:
		case icSigMch7Data:
		case icSigMch8Data: {
				icmLuBase *lu;
				double dval[MAX_CHAN];
				double ncval[3];
				double cvals[MAX_CHAN][3];
				int inn, e, nlighter, ndarker;

				/* Grab a lookup object */
				if ((lu = p->get_luobj(p, icmFwd, icRelativeColorimetric, icSigLabData, icmLuOrdNorm)) == NULL)
					error("icxGetLimits: assert: getting Fwd Lookup failed!");

				lu->spaces(lu, NULL, &inn, NULL, NULL, NULL, NULL, NULL, NULL, NULL);

				/* Decide if the colorspace is aditive or subtractive */

				/* First the no colorant value */
				for (e = 0; e < inn; e++)
					dval[e] = 0.0;
				lu->lookup(lu, ncval, dval);

				/* Then all the colorants */
				nlighter = ndarker = 0;
				for (e = 0; e < inn; e++) {
					dval[e] = 1.0;
					lu->lookup(lu, cvals[e], dval);
					dval[e] = 0.0;
					if (fabs(cvals[e][0] - ncval[0]) > 5.0) {
						if (cvals[e][0] > ncval[0])
							nlighter++;
						else
							ndarker++;
					}
				}

				if (ndarker > 0 && nlighter == 0) {		/* Assume subtractive. */
					double pbk[3] = { 0.0,0.0,0.0 };	/* Perfect black */
					double smd = 1e10;			/* Smallest distance */

					/* Guess the black channel */
					for (e = 0; e < inn; e++) {
						double tt;
						tt = icmNorm33sq(pbk, cvals[e]);
						if (tt < smd) {
							smd = tt;	
							kch = e;
						}
					}
					/* See if the black seems sane */
					if (cvals[kch][0] > 40.0
					 || fabs(cvals[kch][1]) > 10.0
					 || fabs(cvals[kch][2]) > 10.0) {
						kch = -1;
					}
				}
				lu->del(lu);
			}
			break;

		default:
			break;
	}

	return kch;
}

/* Utility function - given an open icc profile, */
/* estmate the total ink limit and black ink limit. */
/* Note that this is rather rough, because ICC profiles */
/* don't have a tag for this information, and ICC profiles */
/* don't have any straightforward way of identifying particular */
/* color channels for > 4 color. */
/* If there are no limits, or they are not discoverable or */
/* applicable, return values of -1.0 */

void icxGetLimits(
xicc *xp,
double *tlimit,
double *klimit 
) {
	icc *p = xp->pp;
	int nch;
	double max[MAX_CHAN];		/* Max of each channel */
	double total;

	total = p->get_tac(p, max, xp->cal != NULL ? xiccCalCallback : NULL, (void *)xp->cal);

	if (total < 0.0) {	/* Not valid */
		if (tlimit != NULL)
			*tlimit = -1.0;
		if (klimit != NULL)
			*klimit = -1.0;
		return;
	}

	nch = icmCSSig2nchan(p->header->colorSpace);

	/* No effective limit */
	if (tlimit != NULL) {
		if (total >= (double)nch) {
			*tlimit = -1.0;
		} else {
			*tlimit = total;
		}
	}

	if (klimit != NULL) {
		int kch;

		kch = icxGuessBlackChan(p);

		if (kch < 0 || max[kch] >= 1.0) {
			*klimit = -1.0;
		} else {
			*klimit = max[kch];
		}
	}
}

/* Replace a non-set limit (ie. < 0.0) with the heuristic from */
/* the given profile. */
void icxDefaultLimits(
xicc *xp,
double *tlout,
double tlin,
double *klout,
double klin
) {
	if (tlin < 0.0 || klin < 0.0) { 
		double tl, kl;

		icxGetLimits(xp, &tl, &kl);

		if (tlin < 0.0)
			tlin = tl;

		if (klin < 0.0)
			klin = kl;
	}

	if (tlout != NULL)
		*tlout = tlin;

	if (klout != NULL)
		*klout = klin;
}

/* Structure to hold optimisation information */
typedef struct {
	xcal *cal;
	double ilimit;
	double uilimit;
} ulimctx;

/* Callback to find equivalent underlying total limit */
/* and try and maximize it while remaining within gamut */
static double ulimitfunc(void *cntx, double pv[]) {
	ulimctx *cx = (ulimctx *)cntx;
	xcal *cal = cx->cal;
	int devchan = cal->devchan;
	int i;
	double dv, odv;
	double og = 0.0, rv = 0.0;

	double usum = 0.0, sum = 0.0;

	/* Comute calibrated sum of channels except last */
	for (i = 0; i < (devchan-1); i++) {
		double dv = pv[i];		/* Underlying (pre-calibration) device value */
		usum += dv;				/* Underlying sum */
		if (dv < 0.0) {
			og += -dv;
			dv = 0.0;
		} else if (dv > 1.0) {
			og += dv - 1.0;
			dv = 1.0;
		} else 
			dv = cal->interp_ch(cal, i, dv);	/* Calibrated device value */
		sum += dv;							/* Calibrated device sum */
	}
	/* Compute the omitted channel value */
	dv = cx->ilimit - sum;					/* Omitted calibrated device value */
	if (dv < 0.0) {
		og += -dv;
		dv = 0.0;
	} else if (dv > 1.0) {
		og += dv - 1.0;
		dv = 1.0;
	} else 
		dv = cal->inv_interp_ch(cal, i, dv);	/* Omitted underlying device value */
	usum += dv;								/* Underlying sum */
	cx->uilimit = usum; 
	
	rv = 10000.0 * og - usum;		/* Penalize out of gamut, maximize underlying sum */

//printf("~1 returning %f from %f %f %f %f\n",rv,pv[0],pv[1],pv[2],dv);
	return rv;
}

/* Given a calibrated total ink limit and an xcal, return the */
/* equivalent underlying (pre-calibration) total ink limit. */
/* This is the maximum equivalent, that makes sure that */
/* the calibrated limit is met or exceeded. */
double icxMaxUnderlyingLimit(xcal *cal, double ilimit) {
	ulimctx cx;
	int i;
	double dv[MAX_CHAN];
	double sr[MAX_CHAN];
	double rv;				/* Residual value */

	if (cal->devchan <= 1) {
		return cal->inv_interp_ch(cal, 0, ilimit);
	}

	cx.cal = cal;
	cx.ilimit = ilimit;

	for (i = 0; i < (cal->devchan-1); i++) {
		sr[i] = 0.05;
		dv[i] = 0.1;
	}
	if (powell(&rv, cal->devchan-1, dv, sr, 0.000001, 1000, ulimitfunc,
			                      (void *)&cx, NULL, NULL) != 0) {
		warning("icxUnderlyingLimit() failed for chan %d, ilimit %f\n",cal->devchan,ilimit);
		return ilimit;
	}
	ulimitfunc((void *)&cx, dv);

	return cx.uilimit;
}

/* ------------------------------------------------------ */
/* Conversion and deltaE formular that include partial */
/* derivatives, for use within fit parameter optimisations. */

/* CIE XYZ to perceptual Lab with partial derivatives. */
void icxdXYZ2Lab(icmXYZNumber *w, double *out, double dout[3][3], double *in) {
	double wp[3], tin[3], dtin[3];
	int i;

	wp[0] = w->X, wp[1] = w->Y, wp[2] = w->Z;

	for (i = 0; i < 3; i++) {
		tin[i] = in[i]/wp[i];
		dtin[i] = 1.0/wp[i];

		if (tin[i] > 0.008856451586) {
			dtin[i] *= pow(tin[i], -2.0/3.0) / 3.0;
			tin[i] = pow(tin[i],1.0/3.0);
		} else {
			dtin[i] *= 7.787036979;
			tin[i] = 7.787036979 * tin[i] + 16.0/116.0;
		}
	}

	out[0] = 116.0 * tin[1] - 16.0;
	dout[0][0] = 0.0;
	dout[0][1] = 116.0 * dtin[1];
	dout[0][2] = 0.0;

	out[1] = 500.0 * (tin[0] - tin[1]);
	dout[1][0] = 500.0 * dtin[0];
	dout[1][1] = 500.0 * -dtin[1];
	dout[1][2] = 0.0;

	out[2] = 200.0 * (tin[1] - tin[2]);
	dout[2][0] = 0.0;
	dout[2][1] = 200.0 * dtin[1];
	dout[2][2] = 200.0 * -dtin[2];
}


/* Return the normal Delta E squared, given two Lab values, */
/* including partial derivatives. */
double icxdLabDEsq(double dout[2][3], double *Lab0, double *Lab1) {
	double rv = 0.0, tt;

	tt = Lab0[0] - Lab1[0];
	dout[0][0] =  2.0 * tt;
	dout[1][0] = -2.0 * tt;
	rv += tt * tt;
	tt = Lab0[1] - Lab1[1];
	dout[0][1] =  2.0 * tt;
	dout[1][1] = -2.0 * tt;
	rv += tt * tt;
	tt = Lab0[2] - Lab1[2];
	dout[0][2] =  2.0 * tt;
	dout[1][2] = -2.0 * tt;
	rv += tt * tt;
	return rv;
}

/* Return the CIE94 Delta E color difference measure, squared */
/* including partial derivatives. */
double icxdCIE94sq(double dout[2][3], double Lab0[3], double Lab1[3]) {
	double desq, _desq[2][3];
	double dlsq;
	double dcsq, _dcsq[2][2];	/* == [x][1,2] */
	double c12,   _c12[2][2];	/* == [x][1,2] */
	double dhsq, _dhsq[2][2];	/* == [x][1,2] */
	double rv;

	{
		double dl, da, db;

		dl = Lab0[0] - Lab1[0];
		dlsq = dl * dl;		/* dl squared */
		da = Lab0[1] - Lab1[1];
		db = Lab0[2] - Lab1[2];


		/* Compute normal Lab delta E squared */
		desq = dlsq + da * da + db * db;
		_desq[0][0] =  2.0 * dl;
		_desq[1][0] = -2.0 * dl;
		_desq[0][1] =  2.0 * da;
		_desq[1][1] = -2.0 * da;
		_desq[0][2] =  2.0 * db;
		_desq[1][2] = -2.0 * db;
	}

	{
		double c1, c2, dc, tt;

		/* Compute chromanance for the two colors */
		c1 = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]);
		c2 = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]);
		c12 = sqrt(c1 * c2);	/* Symetric chromanance */

		tt = 0.5 * (pow(c2, 0.5) + 1e-12)/(pow(c1, 1.5) + 1e-12);
		_c12[0][0] = Lab0[1] * tt;
		_c12[0][1] = Lab0[2] * tt;
		tt = 0.5 * (pow(c1, 0.5) + 1e-12)/(pow(c2, 1.5) + 1e-12);
		_c12[1][0] = Lab1[1] * tt;
		_c12[1][1] = Lab1[2] * tt;

		/* delta chromanance squared */
		dc = c2 - c1;
		dcsq = dc * dc;
		if (c1 < 1e-12 || c2 < 1e-12) {
			c1 += 1e-12;
			c2 += 1e-12;
		}
		_dcsq[0][0] = -2.0 * Lab0[1] * (c2 - c1)/c1;
		_dcsq[0][1] = -2.0 * Lab0[2] * (c2 - c1)/c1;
		_dcsq[1][0] =  2.0 * Lab1[1] * (c2 - c1)/c2;
		_dcsq[1][1] =  2.0 * Lab1[2] * (c2 - c1)/c2;
	}

	/* Compute delta hue squared */
	dhsq = desq - dlsq - dcsq;
	if (dhsq >= 0.0) {
		_dhsq[0][0] = _desq[0][1] - _dcsq[0][0];
		_dhsq[0][1] = _desq[0][2] - _dcsq[0][1];
		_dhsq[1][0] = _desq[1][1] - _dcsq[1][0];
		_dhsq[1][1] = _desq[1][2] - _dcsq[1][1];
	} else {
		dhsq = 0.0;
		_dhsq[0][0] = 0.0;
		_dhsq[0][1] = 0.0;
		_dhsq[1][0] = 0.0;
		_dhsq[1][1] = 0.0;
	}

	{
		double sc, scsq, scf;
		double sh, shsq, shf;

		/* Weighting factors for delta chromanance & delta hue */
		sc = 1.0 + 0.048 * c12;
		scsq = sc * sc;
		
		sh = 1.0 + 0.014 * c12;
		shsq = sh * sh;
	
		rv = dlsq + dcsq/scsq + dhsq/shsq;

		scf = 0.048 * -2.0 * dcsq/(scsq * sc);
		shf = 0.014 * -2.0 * dhsq/(shsq * sh);
		dout[0][0] = _desq[0][0];
		dout[0][1] = _dcsq[0][0]/scsq + _c12[0][0] * scf
		           + _dhsq[0][0]/shsq + _c12[0][0] * shf;
		dout[0][2] = _dcsq[0][1]/scsq + _c12[0][1] * scf
		           + _dhsq[0][1]/shsq + _c12[0][1] * shf;
		dout[1][0] = _desq[1][0];
		dout[1][1] = _dcsq[1][0]/scsq + _c12[1][0] * scf
		           + _dhsq[1][0]/shsq + _c12[1][0] * shf;
		dout[1][2] = _dcsq[1][1]/scsq + _c12[1][1] * scf
		           + _dhsq[1][1]/shsq + _c12[1][1] * shf;
		return rv;
	}
}

// ~~99 not sure if these are correct:

/* Return the normal Delta E given two Lab values, */
/* including partial derivatives. */
double icxdLabDE(double dout[2][3], double *Lab0, double *Lab1) {
	double rv = 0.0, tt;

	tt = Lab0[0] - Lab1[0];
	dout[0][0] =  1.0 * tt;
	dout[1][0] = -1.0 * tt;
	rv += tt * tt;
	tt = Lab0[1] - Lab1[1];
	dout[0][1] =  1.0 * tt;
	dout[1][1] = -1.0 * tt;
	rv += tt * tt;
	tt = Lab0[2] - Lab1[2];
	dout[0][2] =  1.0 * tt;
	dout[1][2] = -1.0 * tt;
	rv += tt * tt;
	return sqrt(rv);
}

/* Return the CIE94 Delta E color difference measure */
/* including partial derivatives. */
double icxdCIE94(double dout[2][3], double Lab0[3], double Lab1[3]) {
	double desq, _desq[2][3];
	double dlsq;
	double dcsq, _dcsq[2][2];	/* == [x][1,2] */
	double c12,   _c12[2][2];	/* == [x][1,2] */
	double dhsq, _dhsq[2][2];	/* == [x][1,2] */
	double rv;

	{
		double dl, da, db;

		dl = Lab0[0] - Lab1[0];
		dlsq = dl * dl;		/* dl squared */
		da = Lab0[1] - Lab1[1];
		db = Lab0[2] - Lab1[2];


		/* Compute normal Lab delta E squared */
		desq = dlsq + da * da + db * db;
		_desq[0][0] =  1.0 * dl;
		_desq[1][0] = -1.0 * dl;
		_desq[0][1] =  1.0 * da;
		_desq[1][1] = -1.0 * da;
		_desq[0][2] =  1.0 * db;
		_desq[1][2] = -1.0 * db;
	}

	{
		double c1, c2, dc, tt;

		/* Compute chromanance for the two colors */
		c1 = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]);
		c2 = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]);
		c12 = sqrt(c1 * c2);	/* Symetric chromanance */

		tt = 0.5 * (pow(c2, 0.5) + 1e-12)/(pow(c1, 1.5) + 1e-12);
		_c12[0][0] = Lab0[1] * tt;
		_c12[0][1] = Lab0[2] * tt;
		tt = 0.5 * (pow(c1, 0.5) + 1e-12)/(pow(c2, 1.5) + 1e-12);
		_c12[1][0] = Lab1[1] * tt;
		_c12[1][1] = Lab1[2] * tt;

		/* delta chromanance squared */
		dc = c2 - c1;
		dcsq = dc * dc;
		if (c1 < 1e-12 || c2 < 1e-12) {
			c1 += 1e-12;
			c2 += 1e-12;
		}
		_dcsq[0][0] = -1.0 * Lab0[1] * (c2 - c1)/c1;
		_dcsq[0][1] = -1.0 * Lab0[2] * (c2 - c1)/c1;
		_dcsq[1][0] =  1.0 * Lab1[1] * (c2 - c1)/c2;
		_dcsq[1][1] =  1.0 * Lab1[2] * (c2 - c1)/c2;
	}

	/* Compute delta hue squared */
	dhsq = desq - dlsq - dcsq;
	if (dhsq >= 0.0) {
		_dhsq[0][0] = _desq[0][1] - _dcsq[0][0];
		_dhsq[0][1] = _desq[0][2] - _dcsq[0][1];
		_dhsq[1][0] = _desq[1][1] - _dcsq[1][0];
		_dhsq[1][1] = _desq[1][2] - _dcsq[1][1];
	} else {
		dhsq = 0.0;
		_dhsq[0][0] = 0.0;
		_dhsq[0][1] = 0.0;
		_dhsq[1][0] = 0.0;
		_dhsq[1][1] = 0.0;
	}

	{
		double sc, scsq, scf;
		double sh, shsq, shf;

		/* Weighting factors for delta chromanance & delta hue */
		sc = 1.0 + 0.048 * c12;
		scsq = sc * sc;
		
		sh = 1.0 + 0.014 * c12;
		shsq = sh * sh;
	
		rv = dlsq + dcsq/scsq + dhsq/shsq;

		scf = 0.048 * -1.0 * dcsq/(scsq * sc);
		shf = 0.014 * -1.0 * dhsq/(shsq * sh);
		dout[0][0] = _desq[0][0];
		dout[0][1] = _dcsq[0][0]/scsq + _c12[0][0] * scf
		           + _dhsq[0][0]/shsq + _c12[0][0] * shf;
		dout[0][2] = _dcsq[0][1]/scsq + _c12[0][1] * scf
		           + _dhsq[0][1]/shsq + _c12[0][1] * shf;
		dout[1][0] = _desq[1][0];
		dout[1][1] = _dcsq[1][0]/scsq + _c12[1][0] * scf
		           + _dhsq[1][0]/shsq + _c12[1][0] * shf;
		dout[1][2] = _dcsq[1][1]/scsq + _c12[1][1] * scf
		           + _dhsq[1][1]/shsq + _c12[1][1] * shf;
		return sqrt(rv);
	}
}

/* ------------------------------------------------------ */
/* A power-like function, based on Graphics Gems adjustment curve. */
/* Avoids "toe" problem of pure power. */
/* Adjusted so that "power" 2 and 0.5 agree with real power at 0.5 */

double icx_powlike(double vv, double pp) {
	double tt, g;

	if (pp >= 1.0) {
		g = 2.0 * (pp - 1.0);
		vv = vv/(g - g * vv + 1.0);
	} else {
		g = 2.0 - 2.0/pp;
		vv = (vv - g * vv)/(1.0 - g * vv);
	}

	return vv;
}

/* Compute the necessary aproximate power, to transform */
/* the given value from src to dst. They are assumed to be */
/* in the range 0.0 .. 1.0 */
double icx_powlike_needed(double src, double dst) {
	double pp, g;

	if (dst <= src) {
		g = -((src - dst)/(dst * src - dst));
		pp = (0.5 * g) + 1.0;
	} else {
		g = -((src - dst)/((dst - 1.0) * src));
		pp = 1.0/(1.0 - 0.5 * g);
	}

	return pp;
}

/* ------------------------------------------------------ */
/* Parameterized transfer/dot gain function. */
/* Used for device modelling. Including partial */
/* derivative for input and parameters. */

/* NOTE that clamping the input values seems to cause */
/* conjgrad() problems. */

/* Transfer function */
double icxTransFunc(
double *v,			/* Pointer to first parameter */
int    luord,		/* Number of parameters */
double vv			/* Source of value */
) {
	double g;
	int ord;

	/* Process all the shaper orders from low to high. */
	/* [These shapers were inspired by a Graphics Gem idea */
	/* (Gems IV, VI.3, "Fast Alternatives to Perlin's Bias and */
	/*  Gain Functions, pp 401). */
	/*  They have the nice properties that they are smooth, and */
	/*  are monotonic. The control parameter has been */
	/*  altered to have a range from -oo to +oo rather than 0.0 to 1.0 */
	/*  so that the search space is less non-linear. */
	for (ord = 0; ord < luord; ord++) {
		int nsec;			/* Number of sections */
		double sec;			/* Section */

		g = v[ord];			/* Parameter */

		nsec = ord + 1;		/* Increase sections for each order */

		vv *= (double)nsec;

		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 /= (double)nsec;
	}

	return vv;
}

/* Inverse transfer function */
double icxInvTransFunc(
double *v,			/* Pointer to first parameter */
int    luord,		/* Number of parameters */
double vv			/* Source of value */
) {
	double g;
	int ord;

	/* Process the shaper orders in reverse from high to low. */
	/* [These shapers were inspired by a Graphics Gem idea */
	/* (Gems IV, VI.3, "Fast Alternatives to Perlin's Bias and */
	/*  Gain Functions, pp 401). */
	/*  They have the nice properties that they are smooth, and */
	/*  are monotonic. The control parameter has been */
	/*  altered to have a range from -oo to +oo rather than 0.0 to 1.0 */
	/*  so that the search space is less non-linear. */
	for (ord = luord-1; ord >= 0; ord--) {
		int nsec;			/* Number of sections */
		double sec;			/* Section */

		g = -v[ord];		/* Inverse parameter */

		nsec = ord + 1;		/* Increase sections for each order */

		vv *= (double)nsec;

		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 /= (double)nsec;
	}

	return vv;
}

/* Transfer function with offset and scale */
double icxSTransFunc(
double *v,			/* Pointer to first parameter */
int    luord,		/* Number of parameters */
double vv,			/* Source of value */
double min,			/* Scale values */
double max
) {
	max -= min;

	vv = (vv - min)/max;
	vv = icxTransFunc(v, luord, vv);
	vv = (vv * max) + min;
	return vv;
}

/* Inverse Transfer function with offset and scale */
double icxInvSTransFunc(
double *v,			/* Pointer to first parameter */
int    luord,		/* Number of parameters */
double vv,			/* Source of value */
double min,			/* Scale values */
double max
) {
	max -= min;

	vv = (vv - min)/max;
	vv = icxInvTransFunc(v, luord, vv);
	vv = (vv * max) + min;
	return vv;
}

/* Transfer function with partial derivative */
/* with respect to the parameters. */
double icxdpTransFunc(
double *v,			/* Pointer to first parameter */
double *dv,			/* Return derivative wrt each parameter */
int    luord,		/* Number of parameters */
double vv			/* Source of value */
) {
	double g;
	int i, ord;

	/* Process all the shaper orders from high to low. */
	for (ord = 0; ord < luord; ord++) {
		double dsv;		/* del for del in g */
		double ddv;		/* del for del in vv */
		int nsec;		/* Number of sections */
		double sec;		/* Section */

		g = v[ord];			/* Parameter */

		nsec = ord + 1;		/* Increase sections for each order */

		vv *= (double)nsec;

		sec = floor(vv);
		if (((int)sec) & 1) {
			g = -g;				/* Alternate action in each section */
		}
		vv -= sec;
		if (g >= 0.0) {
			double tt = g - g * vv + 1.0;
			dsv = (vv * vv - vv)/(tt * tt);
			ddv = (g + 1.0)/(tt * tt);
			vv = vv/tt;
		} else {
			double tt = 1.0 - g * vv;
			dsv = (vv * vv - vv)/(tt * tt);
			ddv = (1.0 - g)/(tt * tt);
			vv = (vv - g * vv)/tt;
		}

		vv += sec;
		vv /= (double)nsec;
		dsv /= (double)nsec;
		if (((int)sec) & 1)
			dsv = -dsv;

		dv[ord] = dsv;
		for (i = ord - 1; i >= 0; i--)
			dv[i] *= ddv;
	}

	return vv;
}

/* Transfer function with offset and scale, and */
/* partial derivative with respect to the parameters. */
double icxdpSTransFunc(
double *v,			/* Pointer to first parameter */
double *dv,			/* Return derivative wrt each parameter */
int    luord,		/* Number of parameters */
double vv,			/* Source of value */
double min,			/* Scale values */
double max
) {
	int i;
	max -= min;

	vv = (vv - min)/max;
	vv = icxdpTransFunc(v, dv, luord, vv);
	vv = (vv * max) + min;
	for (i = 0; i < luord; i++)
		dv[i] *= max;
	return vv;
}


/* Transfer function with partial derivative */
/* with respect to the input value. */
double icxdiTransFunc(
double *v,			/* Pointer to first parameter */
double *pdin,		/* Return derivative wrt source value */
int    luord,		/* Number of parameters */
double vv			/* Source of value */
) {
	double g, din;
	int ord;

#ifdef NEVER
	if (vv < 0.0 || vv > 1.0) {
		if (vv < 0.0)
			vv = 0.0;
	 	else
			vv = 1.0;

		*pdin = 0.0;
		return vv;
	}
#endif
	din = 1.0;

	/* Process all the shaper orders from high to low. */
	for (ord = 0; ord < luord; ord++) {
		double ddv;		/* del for del in vv */
		int nsec;		/* Number of sections */
		double sec;		/* Section */

		g = v[ord];			/* Parameter */

		nsec = ord + 1;		/* Increase sections for each order */

		vv *= (double)nsec;

		sec = floor(vv);
		if (((int)sec) & 1) {
			g = -g;				/* Alternate action in each section */
		}
		vv -= sec;
		if (g >= 0.0) {
			double tt = g - g * vv + 1.0;
			ddv = (g + 1.0)/(tt * tt);
			vv = vv/tt;
		} else {
			double tt = 1.0 - g * vv;
			ddv = (1.0 - g)/(tt * tt);
			vv = (vv - g * vv)/tt;
		}

		vv += sec;
		vv /= (double)nsec;
		din *= ddv;
	}

	*pdin = din;
	return vv;
}

/* Transfer function with offset and scale, and */
/* partial derivative with respect to the input value. */
double icxdiSTransFunc(
double *v,			/* Pointer to first parameter */
double *pdv,		/* Return derivative wrt source value */
int    luord,		/* Number of parameters */
double vv,			/* Source of value */
double min,			/* Scale values */
double max
) {
	max -= min;

	vv = (vv - min)/max;
	vv = icxdiTransFunc(v, pdv, luord, vv);
	vv = (vv * max) + min;
	return vv;
}


/* Transfer function with partial derivative */
/* with respect to the parameters and the input value. */
double icxdpdiTransFunc(
double *v,			/* Pointer to first parameter */
double *dv,			/* Return derivative wrt each parameter */
double *pdin,		/* Return derivative wrt source value */
int    luord,		/* Number of parameters */
double vv			/* Source of value */
) {
	double g, din;
	int i, ord;

#ifdef NEVER
	if (vv < 0.0 || vv > 1.0) {
		if (vv < 0.0)
			vv = 0.0;
	 	else
			vv = 1.0;

		for (ord = 0; ord < luord; ord++)
			dv[ord] = 0.0;
		*pdin = 0.0;
		return vv;
	}
#endif
	din = 1.0;

	/* Process all the shaper orders from high to low. */
	for (ord = 0; ord < luord; ord++) {
		double dsv;		/* del for del in g */
		double ddv;		/* del for del in vv */
		int nsec;		/* Number of sections */
		double sec;		/* Section */

		g = v[ord];			/* Parameter */

		nsec = ord + 1;		/* Increase sections for each order */

		vv *= (double)nsec;

		sec = floor(vv);
		if (((int)sec) & 1) {
			g = -g;				/* Alternate action in each section */
		}
		vv -= sec;
		if (g >= 0.0) {
			double tt = g - g * vv + 1.0;
			dsv = (vv * vv - vv)/(tt * tt);
			ddv = (g + 1.0)/(tt * tt);
			vv = vv/tt;
		} else {
			double tt = 1.0 - g * vv;
			dsv = (vv * vv - vv)/(tt * tt);
			ddv = (1.0 - g)/(tt * tt);
			vv = (vv - g * vv)/tt;
		}

		vv += sec;
		vv /= (double)nsec;
		dsv /= (double)nsec;
		if (((int)sec) & 1)
			dsv = -dsv;

		dv[ord] = dsv;
		for (i = ord - 1; i >= 0; i--)
			dv[i] *= ddv;
		din *= ddv;
	}

	*pdin = din;
	return vv;
}

/* Transfer function with offset and scale, and */
/* partial derivative with respect to the */
/* parameters and the input value. */
double icxdpdiSTransFunc(
double *v,			/* Pointer to first parameter */
double *dv,			/* Return derivative wrt each parameter */
double *pdin,		/* Return derivative wrt source value */
int    luord,		/* Number of parameters */
double vv,			/* Source of value */
double min,			/* Scale values */
double max
) {
	int i;
	max -= min;

	vv = (vv - min)/max;
	vv = icxdpdiTransFunc(v, dv, pdin, luord, vv);
	vv = (vv * max) + min;
	for (i = 0; i < luord; i++)
		dv[i] *= max;
	return vv;
}

/* ------------------------------------------------------ */
/* Multi-plane interpolation, used for device modelling. */
/* Including partial derivative for input and parameters. */
/* A simple flat plane is used for each output. */

/* Multi-plane interpolation - uses base + di slope values.  */
/* Parameters are assumed to be fdi groups of di + 1 parameters. */
void icxPlaneInterp(
double *v,			/* Pointer to first parameter [fdi * (di + 1)] */
int    fdi,			/* Number of output channels */
int    di,			/* Number of input channels */
double *out,		/* Resulting fdi values */
double *in			/* Input di values */
) {
	int e, f;

	for (f = 0; f < fdi; f++) {
		for (out[f] = 0.0, e = 0; e < di; e++, v++) {
			out[f] += in[e] * *v;
		}
		out[f] += *v;
	}
}


/* Multii-plane interpolation with partial derivative */
/* with respect to the input and parameters. */
void icxdpdiPlaneInterp(
double *v,			/* Pointer to first parameter value [fdi * (di + 1)] */
double *dv,			/* Return [1 + di] deriv. wrt each parameter v */
double *din,		/* Return [fdi * di] deriv. wrt each input value */
int    fdi,			/* Number of output channels */
int    di,			/* Number of input channels */
double *out,		/* Resulting fdi values */
double *in			/* Input di values */
) {
	int e, ee, f, g;
	int dip2 = (di + 1);	/* Output dim increment through parameters */

	/* Compute the output values */
	for (f = 0; f < fdi; f++) {
		for (out[f] = 0.0, e = 0; e < di; e++)
			out[f] += in[e] * v[f * dip2 + e];
		out[f] += v[f * dip2 + e];
	}

	/* Since interpolation is verys simple, derivative are also simple */

	/* Copy del for parameter to return array */
	for (e = 0; e < di; e++)
		dv[e] = in[e];
	dv[e] = 1.0;

	/* Compute del of out[] from in[] */
	for (f = 0; f < fdi; f++) {
		for (e = 0; e < di; e++) {
			din[f * di + e] = v[f * dip2 + e];
		}
	}
}

/* ------------------------------------------------------ */
/* Matrix cube interpolation, used for device modelling. */
/* Including partial derivative for input and parameters. */

/* Matrix cube interpolation - interpolate between 2^di output corner values. */
/* Parameters are assumed to be fdi groups of 2^di parameters. */
void icxCubeInterp(
double *v,			/* Pointer to first parameter */
int    fdi,			/* Number of output channels */
int    di,			/* Number of input channels */
double *out,		/* Resulting fdi values */
double *in			/* Input di values */
) {
	int e, f, g;
	double gw[1 << MXDI];		/* weight for each matrix grid cube corner */

	/* Compute corner weights needed for interpolation */
	gw[0] = 1.0;
	for (e = 0, g = 1; e < di; e++, g *= 2) {
		int i; 
		for (i = 0; i < g; i++) {
			gw[g+i] = gw[i] * in[e];
			gw[i] *= (1.0 - in[e]);
		}
	}

	/* Now compute the output values */
	for (f = 0; f < fdi; f++) {
		out[f] = 0.0;						/* For each output value */
		for (e = 0; e < (1 << di); e++) {	/* For all corners of cube */
			out[f] += gw[e] * *v;
			v++;
		}
	}
}


/* Matrix cube interpolation. with partial derivative */
/* with respect to the input and parameters. */
void icxdpdiCubeInterp(
double *v,			/* Pointer to first parameter value [fdi * 2^di] */
double *dv,			/* Return [2^di] deriv. wrt each parameter v */
double *din,		/* Return [fdi * di] deriv. wrt each input value */
int    fdi,			/* Number of output channels */
int    di,			/* Number of input channels */
double *out,		/* Resulting fdi values */
double *in			/* Input di values */
) {
	int e, ee, f, g;
	int dip2 = (1 << di);
	double gw[1 << MXDI];			/* weight for each matrix grid cube corner */

	/* Compute corner weights needed for interpolation */
	gw[0] = 1.0;
	for (e = 0, g = 1; e < di; e++, g *= 2) {
		int i; 
		for (i = 0; i < g; i++) {
			gw[g+i] = gw[i] * in[e];
			gw[i] *= (1.0 - in[e]);
		}
	}

	/* Now compute the output values */
	for (f = 0; f < fdi; f++) {
		out[f] = 0.0;						/* For each output value */
		for (ee = 0; ee < dip2; ee++) {	/* For all corners of cube */
			out[f] += gw[ee] * v[f * dip2 + ee];
		}
	}

	/* Copy del for parameter to return array */
	for (ee = 0; ee < dip2; ee++) {	/* For all other corners of cube */
		dv[ee] = gw[ee];					/* del from parameter */
	}

	/* Compute del from in[] value we want */
	for (e = 0; e < di; e++) {				/* For input we want del wrt */

		for (f = 0; f < fdi; f++)
			din[f * di + e] = 0.0;			/* Zero del ready of accumulation */

		for (ee = 0; ee < dip2; ee++) {		/* For all corners of cube weights, */
			int e2;							/* accumulate del from in[] we want. */
			double vv = 1.0;

			/* Compute in[] weighted cube corners for all except del of in[] we want */
			for (e2 = 0; e2 < di; e2++) {			/* const from non del inputs */
				if (e2 == e)
					continue;
				if (ee & (1 << e2))
					vv *= in[e2];
				else
					vv *= (1.0 - in[e2]);
			}

			/* Accumulate contribution of in[] we want for corner to out[] we want */
			if (ee & (1 << e)) {
				for (f = 0; f < fdi; f++)
					din[f * di + e] += v[f * dip2 + ee] * vv;
			} else {
				for (f = 0; f < fdi; f++)
					din[f * di + e] -= v[f * dip2 + ee] * vv;
			}
		}
	}
}

/* ------------------------------------------------------ */
/* Matrix cube simplex interpolation, used for device modelling. */

/* Matrix cube simplex interpolation - interpolate between 2^di output corner values. */
/* Parameters are assumed to be fdi groups of 2^di parameters. */
void icxCubeSxInterp(
double *v,			/* Pointer to first parameter */
int    fdi,			/* Number of output channels */
int    di,			/* Number of input channels */
double *out,		/* Resulting fdi values */
double *in			/* Input di values */
) {
	int    si[MAX_CHAN];		/* in[] Sort index, [0] = smalest */

//{
//	double tout[MXDO];
//
//	icxCubeInterp(v, fdi, di, tout, in);
//printf("\n~1 Cube interp result = %f\n",tout[0]);
//}

//printf("~1 icxCubeSxInterp: %f %f %f\n", in[0], in[1], in[2]);
	/* Do insertion sort on coordinates, smallest to largest. */
	{
		int ff, vf;
		unsigned int e;
		double v;
		for (e = 0; e < di; e++)
			si[e] = e;						/* Initial unsorted indexes */

		for (e = 1; e < di; e++) {
			ff = e;
			v = in[si[ff]];
			vf = ff;
			while (ff > 0 && in[si[ff-1]] > v) {
				si[ff] = si[ff-1];
				ff--;
			}
			si[ff] = vf;
		}
	}
//printf("~1 sort order %d %d %d\n", si[0], si[1], si[2]);
//printf("         from %f %f %f\n", in[si[0]], in[si[1]], in[si[2]]);

	/* Now compute the weightings, simplex vertices and output values */
	{
		unsigned int e, f;
		double w;		/* Current vertex weight */

		w = 1.0 - in[si[di-1]];		/* Vertex at base of cell */
		for (f = 0; f < fdi; f++) {
			out[f] = w * v[f * (1 << di)];
//printf("~1 out[%d] = %f = %f * %f\n",f,out[f],w,v[f * (1 << di)]);
		}

		for (e = di-1; e > 0; e--) {	/* Middle verticies */
			w = in[si[e]] - in[si[e-1]];
			v += (1 << si[e]);				/* Move to top of cell in next largest dimension */
			for (f = 0; f < fdi; f++) {
				out[f] += w * v[f * (1 << di)];
//printf("~1 out[%d] = %f += %f * %f\n",f,out[f],w,v[f * (1 << di)]);
			}
		}

		w = in[si[0]];
		v += (1 << si[0]);		/* Far corner from base of cell */
		for (f = 0; f < fdi; f++) {
			out[f] += w * v[f * (1 << di)];
//printf("~1 out[%d] = %f += %f * %f\n",f,out[f],w,v[f * (1 << di)]);
		}
	}
}

/* ------------------------------------------------------ */
/* Matrix multiplication, used for device modelling. */
/* Including partial derivative for input and parameters. */


/* 3x3 matrix multiplication, with the matrix in a 1D array */
/* with respect to the input and parameters. */
void icxMulBy3x3Parm(
	double out[3],			/* Return input multiplied by matrix */
	double mat[9],			/* Matrix organised in [slow][fast] order */
	double in[3]			/* Input values */
) {
	double *v = mat, ov[3];
	int e, f;

	/* Compute the output values */
	for (f = 0; f < 3; f++) {
		ov[f] = 0.0;						/* For each output value */
		for (e = 0; e < 3; e++) {
			ov[f] += *v++ * in[e];
		}
	}
	out[0] = ov[0];
	out[1] = ov[1];
	out[2] = ov[2];
}


/* 3x3 matrix multiplication, with partial derivatives */
/* with respect to the input and parameters. */
void icxdpdiMulBy3x3Parm(
	double out[3],			/* Return input multiplied by matrix */
	double dv[3][9],		/* Return deriv for each [output] with respect to [param] */
	double din[3][3],		/* Return deriv for each [output] with respect to [input] */
	double mat[9],			/* Matrix organised in [slow][fast] order */
	double in[3]			/* Input values */
) {
	double *v, ov[3];
	int e, f;

	/* Compute the output values */
	v = mat;
	for (f = 0; f < 3; f++) {
		ov[f] = 0.0;						/* For each output value */
		for (e = 0; e < 3; e++) {
			ov[f] += *v++ * in[e];
		}
	}

	/* Compute deriv. with respect to the matrix parameter % 3 */
	/* This is pretty simple for a matrix ... */
	for (f = 0; f < 3; f++) {
		for (e = 0; e < 9; e++) {
			if (e/3 ==  f)
				dv[f][e] = in[e % 3];
			else
				dv[f][e] = 0.0;
		}
	}
	
	/* Compute deriv. with respect to the input values */
	/* This is pretty simple for a matrix ... */
	v = mat;
	for (f = 0; f < 3; f++)
		for (e = 0; e < 3; e++)
			din[f][e] = *v++;

	out[0] = ov[0];
	out[1] = ov[1];
	out[2] = ov[2];
}

#undef stricmp