/* 
 * International Color Consortium color transform expanded support
 *
 * Author:  Graeme W. Gill
 * Date:    1/5/01
 * Version: 1.20
 *
 * Copyright 2000-2004 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.
 *
 */

/*
 * This is some test code to test the CIECAM97s3 functionality. 
 */


#include <stdio.h>
#include <math.h>
#include "xcam.h"
#include "cam97s3.h"

#define DIAG		/* Print internal value diagnostics for each conversion */
#define SPOTTEST	/* Test known spot colors */
#undef INVTEST		/* Lab cube to XYZ to Jab to XYZ */
#undef INVTEST2	/* Jab cube to XYZ to Jab */
#define TRES 33
#define USE_HK 0	/* Use Helmholtz-Kohlraush */
#define NOCAMCLIP 1	/* Don't clip before XYZ2Jab */

#ifndef _isnan
#define _isnan(x) ((x) != (x))
#define _finite(x) ((x) == (x))
#endif

static void
Lab2XYZ(double *out, double *in) {
	double L = in[0], a = in[1], b = in[2];
	double x,y,z,fx,fy,fz;

	if (L > 8.0) {
		fy = (L + 16.0)/116.0;
		y = pow(fy,3.0);
	} else {
		y = L/903.2963058;
		fy = 7.787036979 * y + 16.0/116.0;
	}

	fx = a/500.0 + fy;
	if (fx > 24.0/116.0)
		x = pow(fx,3.0);
	else
		x = (fx - 16.0/116.0)/7.787036979;

	fz = fy - b/200.0;
	if (fz > 24.0/116.0)
		z = pow(fz,3.0);
	else
		z = (fz - 16.0/116.0)/7.787036979;

	out[0] = x * 0.9642;
	out[1] = y * 1.0000;
	out[2] = z * 0.8249;
}

/* CIE XYZ to perceptual Lab */
static void
XYZ2Lab(double *out, double *in) {
	double X = in[0], Y = in[1], Z = in[2];
	double x,y,z,fx,fy,fz;

	x = X/0.9642;
	y = Y/1.0000;
	z = Z/0.8249;

	if (x > 0.008856451586)
		fx = pow(x,1.0/3.0);
	else
		fx = 7.787036979 * x + 16.0/116.0;

	if (y > 0.008856451586)
		fy = pow(y,1.0/3.0);
	else
		fy = 7.787036979 * y + 16.0/116.0;

	if (z > 0.008856451586)
		fz = pow(z,1.0/3.0);
	else
		fz = 7.787036979 * z + 16.0/116.0;

	out[0] = 116.0 * fy - 16.0;
	out[1] = 500.0 * (fx - fy);
	out[2] = 200.0 * (fy - fz);
}


/* Return maximum difference */
double maxdiff(double in1[3], double in2[3]) {
	int i;
	double tt, rv = 0.0;

	/* Deal with the possibility that we have nan's */
	for (i = 0; i < 3; i++) {
		tt = fabs(in1[i] - in2[i]);
		if (_isnan(tt))
			return tt;
		if (tt > rv)
			rv = tt;
	}
	return rv;
}

/* Return absolute difference */
double absdiff(double in1[3], double in2[3]) {
	double tt, rv = 0.0;
	tt = in1[0] - in2[0];
	rv += tt * tt;
	tt = in1[1] - in2[1];
	rv += tt * tt;
	tt = in1[2] - in2[2];
	rv += tt * tt;
	return sqrt(rv);
}

/* Return maximum Lab difference of XYZ */
double maxxyzdiff(double i1[3], double i2[3]) {
	int i;
	double tt, rv = 0.0;
	double in1[3], in2[3];

	XYZ2Lab(in1, i1);
	XYZ2Lab(in2, i2);

	/* Deal with the possibility that we have nan's */
	for (i = 0; i < 3; i++) {
		tt = fabs(in1[i] - in2[i]);
		if (_isnan(tt))
			return tt;
		if (tt > rv)
			rv = tt;
	}
	return rv;
}

int
main(void) {
	int ok = 1;
	double white[6][3] = {
		{ 0.9505, 1.000, 1.0888 },
		{ 0.9505, 1.000, 1.0888 },
		{ 1.0985, 1.000, 0.3558 },
		{ 1.0985, 1.000, 0.3558 },
		{ 0.9505, 1.000, 1.0888 },
		{ 0.9642, 1.000, 0.8249 }	/* D50 for inversion tests */
	};
	double La[6] = { 318.31, 31.83, 318.31, 31.83, 318.31, 150.0 };
	double sample[5][3] = {
		{ 0.1901, 0.2000, 0.2178 },
		{ 0.5706, 0.4306, 0.3196 },
		{ 0.0353, 0.0656, 0.0214 },
		{ 0.1901, 0.2000, 0.2178 },
		{ 0.9505, 1.000, 1.0888 }		/* Check white */
	};

#ifdef CIECAM97S3_HK
#ifdef CIECAM97S3_SPLINE_E
	double correct[5][4] = {			/* Hue range is last two */
		{ 41.14,     0.06,     225.3, 251.9 },
		{ 69.04,     71.04,    19.4,  19.4 },
		{ 35.09,     87.15,   175.3, 175.3 },
		{ 55.64,     82.43,   252.5, 252.5 },
		{ 100.00000, 0.0,     0.0,   360.0 }
	};
#else
	double correct[5][4] = {			/* Hue range is last two */
		{ 41.14,     0.06,     225.3, 251.9 },
		{ 69.056048, 71.20763, 19.4,  19.4 },
		{ 35.365009, 88.65238, 175.3, 175.3 },
		{ 55.266184, 80.55128, 252.5, 252.5 },
		{ 100.00000, 0.0,      0.0,   360.0 }
	};
#endif
#else
	double correct[5][4] = {			/* Hue range is last two */
		{ 41.13,  0.05,   225.3, 251.9 },
		{ 64.14,  71.22,  19.4,  19.4 },
		{ 19.18,  88.64,  175.3, 175.3 },
		{ 39.11,  80.55,  252.5, 252.5 },
		{ 100.0,  0.00,   0.0,   360.0 }		/* Check white */
	};
#endif

	double Jab[3];
	double res[3];
	cam97s3 *cam;
	int c;

	cam = new_cam97s3();

#ifdef SPOTTEST
	for (c = 0; c < 5; c++) {

#ifdef DIAG
		printf("Case %d:\n",c);
#endif /* DIAG */
		cam->set_view(
			cam,
			vc_average,	/* Enumerated Viewing Condition */
			white[c],	/* Reference/Adapted White XYZ (Y range 0.0 .. 1.0) */
			La[c],		/* Adapting/Surround Luminance cd/m^2 */
			0.20,		/* Relative Luminance of Background to reference white */
			0.0,		/* Luminance of white in image - not used */
			0.00,		/* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
			white[c],	/* The Flare color coordinates (typically the Ambient color) */
			USE_HK
		);
#ifdef DIAG
		printf("\n");
#endif /* DIAG */
	
		{
			double JCh[3], target[3];
			cam->XYZ_to_cam(cam, Jab, sample[c]);
	
			/* Convert to JCh for checking */
			JCh[0] = Jab[0];

			/* Compute hue angle */
    		JCh[2] = (180.0/3.14159265359) * atan2(Jab[2], Jab[1]);
			JCh[2] = (JCh[2] < 0.0) ? JCh[2] + 360.0 : JCh[2];
	
			/* Compute chroma value */
			JCh[1] = sqrt(Jab[1] * Jab[1] + Jab[2] * Jab[2]);

			target[0] = correct[c][0];
			target[1] = correct[c][1];
			if (JCh[2] >= correct[c][2] && JCh[2] <= correct[c][3]) {
				target[2] = JCh[2];
			} else {
				if (fabs(JCh[2] - correct[c][2]) < fabs(JCh[2] - correct[c][3]))
					target[2] = correct[c][2];
				else
					target[2] = correct[c][3];
			}

			cam->cam_to_XYZ(cam, res, Jab);

			if (maxdiff(JCh, target) > 0.05) {
				printf("Spottest: Excesive error in conversion to CAM %f\n",maxdiff(JCh, target));
				ok = 0;
			}
			if (maxdiff(sample[c], res) > 1e-5) {
				printf("Spottest: Excessive error in inversion %f\n",maxdiff(sample[c], res));
				ok = 0;
			}
#ifdef DIAG
			printf("Jab is %f %f %f, Jch is %f %f %f\n",
			        Jab[0], Jab[1], Jab[2],
			        JCh[0], JCh[1], JCh[2]);

			printf("Error to expected value = %f\n",maxdiff(JCh, target));
		
			printf("XYZ is %f %f %f\n",res[0], res[1], res[2]);
			printf("Inversion error = %f\n",maxdiff(sample[c], res));
#endif /* DIAG */
		}
	}
#endif /* SPOTTEST */

#ifdef INVTEST
	for (c = 5; c < 6; c++) {
		/* Get the range of Jab values */
		double min[3] = { 1e6, 1e6, 1e6 };
		double max[3] = { -1e6, -1e6, -1e6 };

		cam->set_view(
			cam,
			vc_average,	/* Enumerated Viewing Condition */
			white[c],	/* Reference/Adapted White XYZ (Y range 0.0 .. 1.0) = D50 */
			34.0,		/* Adapting/Surround Luminance cd/m^2 */
			0.20,		/* Relative Luminance of Background to reference white */
			0.0,		/* Luminance of white in image - not used */
			0.01,		/* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
			white[c],	/* The Flare color coordinates (typically the Ambient color) */
			USE_HK
		);
	
		{
			int co0, co1, co2;		/* (using co[3] triggers compiler bug) */
			double merr = 0.0;
			double xyz[3], Lab[3], Jab[3], checkxyz[3];


#ifdef NEVER
			/* Test case */
			Lab[0] = 0.0;
			Lab[1] = -180.0;
			Lab[2] = -180.0;
			Lab2XYZ(xyz, Lab);
			cam->XYZ_to_cam(cam, Jab, xyz);
			cam->cam_to_XYZ(cam, checkxyz, Jab);

#else /* !NEVER */

			/* itterate through Lab space */
			for (co0 = 0; co0 < TRES; co0++) {
				Lab[0] = co0/(TRES-1.0);
				Lab[0] = Lab[0] * 100.0;
				for (co1 = 0; co1 < TRES; co1++) {
					Lab[1] = co1/(TRES-1.0);
					Lab[1] = (Lab[1] - 0.5) * 256.0;
					for (co2 = 0; co2 < TRES; co2++) {
						double mxd;
						int i;
						Lab[2] = co2/(TRES-1.0);
						Lab[2] = (Lab[2] - 0.5) * 256.0;
		
						Lab2XYZ(xyz, Lab);
						cam->XYZ_to_cam(cam, Jab, xyz);
						for (i = 0; i < 3; i++) {
							if (Jab[i] < min[i])
								min[i] = Jab[i];
							if (Jab[i] > max[i])
								max[i] = Jab[i];
						}
						cam->cam_to_XYZ(cam, checkxyz, Jab);

						/* Check the result */
						mxd = maxxyzdiff(checkxyz, xyz);
						if (_finite(merr) && (_isnan(mxd) || mxd > merr))
							merr = mxd;
#ifdef DIAG
						if (_isnan(mxd) || mxd > 0.1) {
							printf("\n");
							printf("Lab = %f %f %f\n",Lab[0], Lab[1], Lab[2]);
							printf("%f %f %f -> %f %f %f -> %f %f %f [%f]\n",
							xyz[0],xyz[1],xyz[2],Jab[0],Jab[1],Jab[2],
							checkxyz[0],checkxyz[1],checkxyz[2],
							maxxyzdiff(checkxyz, xyz));
						}
#endif /* DIAG */
					}
				}
			}
			if (_isnan(merr) || merr > 0.1) {
				printf("Excessive error in inversion check %f\n",merr);
				ok = 0;
			}
#ifdef DIAG
			printf("Inversion check complete, peak error = %f\n",merr);
#endif /* !NEVER */

#endif /* DIAG */
			printf("\n");
			printf("Range of Jab values was:\n");
			printf("J:  %f -> %f\n", min[0], max[0]);
			printf("a:  %f -> %f\n", min[1], max[1]);
			printf("b:  %f -> %f\n", min[2], max[2]);
		}
	}
#endif /* INVTEST */

#ifdef INVTEST2
	for (c = 5; c < 6; c++) {
		cam->set_view(
			cam,
			vc_average,	/* Enumerated Viewing Condition */
			white[c],	/* Reference/Adapted White XYZ (Y range 0.0 .. 1.0) = D50 */
			34.0,		/* Adapting/Surround Luminance cd/m^2 */
			0.20,		/* Relative Luminance of Background to reference white */
			0.0,		/* Luminance of white in image - not used */
			0.01,		/* Flare as a fraction of the reference white (Y range 0.0 .. 1.0) */
			white[c],	/* The Flare color coordinates (typically the Ambient color) */
			USE_HK
		);
	
		{
			int co0, co1, co2;		/* (using co[3] triggers compiler bug) */
			double merr = 0.0;
			double xyz[3], Jab[3], checkJab[3];

#ifdef NEVER
			/* Test case */
#ifdef NEVER
			Jab[0] = 0.0;
			Jab[1] = -180.0;
			Jab[2] = -180.0;
			cam->cam_to_XYZ(cam, xyz, Jab);
			cam->XYZ_to_cam(cam, checkJab, xyz);
#else
			xyz[0] = 555.365115;
			xyz[1] = -518.091908;
			xyz[2] = -1591.731701;
			cam->XYZ_to_cam(cam, Jab, xyz);
			cam->cam_to_XYZ(cam, xyz, Jab);
#endif

#else /* !NEVER */

			/* itterate through Jab space */
			for (co0 = 0; co0 < TRES; co0++) {
				Jab[0] = co0/(TRES-1.0);
				Jab[0] = 0.0 + Jab[0] * 100.0;	
				for (co1 = 0; co1 < TRES; co1++) {
					Jab[1] = co1/(TRES-1.0);
					Jab[1] = (Jab[1] - 0.5) * 256.0;
					for (co2 = 0; co2 < TRES; co2++) {
						double mxd;
						Jab[2] = co2/(TRES-1.0);
						Jab[2] = (Jab[2] - 0.5) * 256.0;
		
						cam->cam_to_XYZ(cam, xyz, Jab);
						cam->XYZ_to_cam(cam, checkJab, xyz);

						/* Check the result */
						mxd = maxdiff(checkJab, Jab);
						if (_finite(merr) && (_isnan(mxd) || mxd > merr))
							merr = mxd;
#ifdef DIAG
						if (_isnan(mxd) || mxd > 0.1) {
							printf("\n");
							printf("Jab = %f %f %f\n",Jab[0], Jab[1], Jab[2]);
							printf("%f %f %f -> %f %f %f -> %f %f %f [%f]\n",
							Jab[0],Jab[1],Jab[2], xyz[0],xyz[1],xyz[2],
							checkJab[0],checkJab[1],checkJab[2],
							maxdiff(checkJab, Jab));
						}
#endif /* DIAG */
					}
				}
			}
			if (_isnan(merr) || merr > 1.0) {
				printf("Excessive error in inversion check 2 %f\n",merr);
				ok = 0;
			}
#ifdef DIAG
			printf("Inversion check 2 complete, peak error = %f\n",merr);
#endif /* DIAG */
#endif /* !NEVER */
		}
	}
#endif /* INVTEST2 */

	printf("\n");
	if (ok == 0) {
		printf("Cam testing FAILED\n");
	} else {
		printf("Cam testing OK\n");
	}

	return 0;
}