/* 
 * Argyll Color Correction System
 * Color Device profile checker.
 *
 * Author: Graeme W. Gill
 * Date:   15/7/2001
 *
 * Copyright 2001 - 2005 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 program takes in the .ti3 scattered test chart
 * points, and checks them against an ICC profile.
 * forward ICC device profile.
 */

/*
 * TTBD:
 *		Switch to generic colorant read code rather than Grey/RGB/CMYK,
 *		and allow checking ICC profiles > 4 colors
 */

#undef DEBUG

#define IMP_MONO			/* Turn on development code */

#define verbo stdout

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include "copyright.h"
#include "aconfig.h"
#include "numlib.h"
#include "cgats.h"
#include "xicc.h"
#include "insttypes.h"
#include "sort.h"

void
usage(void) {
	fprintf(stderr,"Check accuracy of ICC profile, Version %s\n",ARGYLL_VERSION_STR);
	fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
	fprintf(stderr,"usage: profcheck [-options] data.ti3 iccprofile.icm\n");
	fprintf(stderr," -v [level]      Verbosity level (default 1), 2 to print each DE\n");
	fprintf(stderr," -c              Show CIE94 delta E values\n");
	fprintf(stderr," -k              Show CIEDE2000 delta E values\n");
	fprintf(stderr," -w              create VRML visualisation (iccprofile.wrl)\n");
	fprintf(stderr," -x              Use VRML axes\n");
	fprintf(stderr," -m              Make VRML lines a minimum of 0.5\n");
	fprintf(stderr," -e              Color vectors acording to delta E\n");
	fprintf(stderr," -d devval1,deval2,devvalN\n");
	fprintf(stderr,"                 Specify a device value to sort against\n");
	fprintf(stderr," -p              Sort device value by PCS (Lab) target\n");
	fprintf(stderr," -f [illum]      Use Fluorescent Whitening Agent compensation [opt. simulated inst. illum.:\n");
	fprintf(stderr,"                  M0, M1, M2, A, C, D50 (def.), D50M2, D65, F5, F8, F10 or file.sp]\n");
	fprintf(stderr," -i illum        Choose illuminant for computation of CIE XYZ from spectral data & FWA:\n");
	fprintf(stderr,"                  A, C, D50 (def.), D50M2, D65, F5, F8, F10 or file.sp\n");
	fprintf(stderr," -o observ       Choose CIE Observer for spectral data:\n");
	fprintf(stderr,"                  1931_2 (def), 1964_10, S&B 1955_2, shaw, J&V 1978_2\n");
	fprintf(stderr," -I intent       r = relative colorimetric, a = absolute (default)\n");
	fprintf(stderr," data.ti3        Test data file\n");
	fprintf(stderr," iccprofile.icm  Profile to check against\n");
	exit(1);
	}

FILE *start_vrml(char *name, int doaxes);
void start_line_set(FILE *wrl);
void add_vertex(FILE *wrl, double pp[3]);
void make_lines(FILE *wrl, int ppset);
void make_de_lines(FILE *wrl);
void end_vrml(FILE *wrl);

/* Patch value type */
typedef struct {
	char sid[50];		/* sample id */
	char slo[50];		/* sample location, "" if not known */
	double p[MAX_CHAN];	/* Device value */
	double v[3];		/* CIE value */
	double dp;			/* Delta from target value */
	double dv;			/* Delta from CIE value */
} pval;

int main(int argc, char *argv[])
{
	int fa,nfa;				/* current argument we're looking at */
	int verb = 0;
	int cie94 = 0;
	int cie2k = 0;
	int dovrml = 0;
	int dominl = 0;
	int doaxes = 0;
	int dodecol = 0;
	char ti3name[MAXNAMEL+1] = { 0 };	/* Input cgats file base name */
	cgats *icg;				/* input cgats structure */
	char iccname[MAXNAMEL+1] = { 0 };	/* Input icc file base name */
	icmFile *rd_fp;
	icRenderingIntent intent = icAbsoluteColorimetric;
	icc *rd_icco;
	icmLuBase *luo;
	char out_name[MAXNAMEL+1], *xl;		/* VRML name */
	FILE *wrl = NULL;

	int fwacomp = 0;			/* FWA compensation */
	int isdisp = 0;				/* nz if this is a display device, 0 if output */
	int isdnormed = 0;      	/* Has display data been normalised to 100 ? */
	int spec = 0;				/* Use spectral data flag */
	icxIllumeType tillum = icxIT_none;	/* Target/simulated instrument illuminant */ 
	xspect cust_tillum, *tillump = NULL; /* Custom target/simulated illumination spectrum */
	icxIllumeType illum = icxIT_D50;	/* Spectral defaults */
	xspect cust_illum;			/* Custom illumination spectrum */
	icxObserverType observ = icxOT_CIE_1931_2;

	int ddevv = 0;				/* Do device value sort */
	double devval[MAX_CHAN];	/* device value to sort on */
	int sortbypcs = 0;			/* Sort by PCS */

	int npat;					/* Number of patches */
	pval *tpat;					/* Patch input values */
	int i, j, rv = 0;
	icColorSpaceSignature devspace = 0;	/* The device colorspace */
	int isAdditive = 0;			/* 0 if subtractive, 1 if additive colorspace */
	int isLab = 0;				/* 0 if input is XYZ, 1 if input is Lab */
	int devchan = 0;			/* Number of device chanels */

#if defined(__IBMC__)
	_control87(EM_UNDERFLOW, EM_UNDERFLOW);
	_control87(EM_OVERFLOW, EM_OVERFLOW);
#endif

	error_program = "profcheck";

	if (argc <= 1)
		usage();

	/* Process the arguments */
	for(fa = 1;fa < argc;fa++) {

		nfa = fa;					/* skip to nfa if next argument is used */
		if (argv[fa][0] == '-') {	/* Look for any flags */
			char *na = NULL;		/* next argument after flag, null if none */

			if (argv[fa][2] != '\000')
				na = &argv[fa][2];		/* next is directly after flag */
			else {
				if ((fa+1) < argc) {
					if (argv[fa+1][0] != '-') {
						nfa = fa + 1;
						na = argv[nfa];		/* next is seperate non-flag argument */
					}
				}
			}

			if (argv[fa][1] == '?')
				usage();

			/* Verbosity */
			else if (argv[fa][1] == 'v') {
				verb = 1;
				if (na != NULL && isdigit(na[0])) {
					verb = atoi(na);
				}
			}

			/* VRML */
			else if (argv[fa][1] == 'w')
				dovrml = 1;

			/* Minimum line length */
			else if (argv[fa][1] == 'm')
				dominl = 1;

			/* Axes */
			else if (argv[fa][1] == 'x')
				doaxes = 1;

			/* Delta E coloring */
			else if (argv[fa][1] == 'e')
				dodecol = 1;

			else if (argv[fa][1] == 'c') {
				cie94 = 1;
				cie2k = 0;
			}

			else if (argv[fa][1] == 'k') {
				cie94 = 0;
				cie2k = 1;
			}

			/* Device sort value */
			else if (argv[fa][1] == 'd') {
				char *tp, buf[200];
				int ndv;
				fa = nfa;
				if (na == NULL) usage();

				ddevv = 1;
				strcpy(buf, na);
				
				/* Replace ',' with '\000' */
				for (ndv = 1,tp = buf; *tp != '\000'; tp++) {
					if (*tp == ',') {
						*tp = '\000';
						ndv++;
					}
				}
				if (ndv >= MAX_CHAN)
					ndv = MAX_CHAN;

				for (tp = buf, i = 0; i < ndv; i++, tp += strlen(tp) + 1) {
					devval[i] = atof(tp);
				}
			}

			else if (argv[fa][1] == 'p')
				sortbypcs = 1;

			/* FWA compensation */
			else if (argv[fa][1] == 'f') {
				fwacomp = 1;

				if (na != NULL) {	/* Argument is present - target/simulated instr. illum. */
					fa = nfa;
					if (strcmp(na, "A") == 0
					 || strcmp(na, "M0") == 0) {
						spec = 1;
						tillum = icxIT_A;
					} else if (strcmp(na, "C") == 0) {
						spec = 1;
						tillum = icxIT_C;
					} else if (strcmp(na, "D50") == 0
					        || strcmp(na, "M1") == 0) {
						spec = 1;
						tillum = icxIT_D50;
					} else if (strcmp(na, "D50M2") == 0
					        || strcmp(na, "M2") == 0) {
						spec = 1;
						tillum = icxIT_D50M2;
					} else if (strcmp(na, "D65") == 0) {
						spec = 1;
						tillum = icxIT_D65;
					} else if (strcmp(na, "F5") == 0) {
						spec = 1;
						tillum = icxIT_F5;
					} else if (strcmp(na, "F8") == 0) {
						spec = 1;
						tillum = icxIT_F8;
					} else if (strcmp(na, "F10") == 0) {
						spec = 1;
						tillum = icxIT_F10;
					} else {	/* Assume it's a filename */
						spec = 1;
						tillum = icxIT_custom;
						if (read_xspect(&cust_tillum, na) != 0)
							usage();
					}
				}
			}
			/* Spectral Illuminant type */
			else if (argv[fa][1] == 'i') {
				fa = nfa;
				if (na == NULL) usage();
				if (strcmp(na, "A") == 0) {
					spec = 1;
					illum = icxIT_A;
				} else if (strcmp(na, "C") == 0) {
					spec = 1;
					illum = icxIT_C;
				} else if (strcmp(na, "D50") == 0) {
					spec = 1;
					illum = icxIT_D50;
				} else if (strcmp(na, "D50M2") == 0) {
					spec = 1;
					illum = icxIT_D50M2;
				} else if (strcmp(na, "D65") == 0) {
					spec = 1;
					illum = icxIT_D65;
				} else if (strcmp(na, "F5") == 0) {
					spec = 1;
					illum = icxIT_F5;
				} else if (strcmp(na, "F8") == 0) {
					spec = 1;
					illum = icxIT_F8;
				} else if (strcmp(na, "F10") == 0) {
					spec = 1;
					illum = icxIT_F10;
				} else {	/* Assume it's a filename */
					spec = 1;
					illum = icxIT_custom;
					if (read_xspect(&cust_illum, na) != 0)
						usage();
				}
			}

			/* Spectral Observer type */
			else if (argv[fa][1] == 'o') {
				fa = nfa;
				if (na == NULL) usage();
				if (strcmp(na, "1931_2") == 0) {			/* Classic 2 degree */
					spec = 1;
					observ = icxOT_CIE_1931_2;
				} else if (strcmp(na, "1964_10") == 0) {	/* Classic 10 degree */
					spec = 1;
					observ = icxOT_CIE_1964_10;
				} else if (strcmp(na, "1955_2") == 0) {		/* Stiles and Burch 1955 2 degree */
					spec = 1;
					observ = icxOT_Stiles_Burch_2;
				} else if (strcmp(na, "1978_2") == 0) {		/* Judd and Voss 1978 2 degree */
					spec = 1;
					observ = icxOT_Judd_Voss_2;
				} else if (strcmp(na, "shaw") == 0) {		/* Shaw and Fairchilds 1997 2 degree */
					spec = 1;
					observ = icxOT_Shaw_Fairchild_2;
				} else
					usage();
			}

			/* Intent (only applies to ICC profile) */
			else if (argv[fa][1] == 'I') {
				fa = nfa;
				if (na == NULL) usage();
    			switch (na[0]) {
					case 'r':
						intent = icRelativeColorimetric;
						break;
					case 'a':
						intent = icAbsoluteColorimetric;
						break;
					default:
						usage();
				}
			}

			else 
				usage();
		} else
			break;
	}

	/* Get the file name arguments */
	if (fa >= argc || argv[fa][0] == '-') usage();
	strncpy(ti3name,argv[fa++],MAXNAMEL); ti3name[MAXNAMEL] = '\000';

	if (fa >= argc || argv[fa][0] == '-') usage();
	strncpy(iccname,argv[fa++],MAXNAMEL); iccname[MAXNAMEL] = '\000';

	strncpy(out_name,iccname,MAXNAMEL-4); out_name[MAXNAMEL-4] = '\000';
	if ((xl = strrchr(out_name, '.')) == NULL)	/* Figure where extention is */
		xl = out_name + strlen(out_name);
	strcpy(xl,".wrl");

	if (fwacomp && spec == 0)
		error ("FWA compensation only works when viewer and/or illuminant selected");

	/* Open and look at the .ti3 profile patches file */
	icg = new_cgats();			/* Create a CGATS structure */
	icg->add_other(icg, "CTI3"); 	/* our special input type is Calibration Target Information 3 */

	if (icg->read_name(icg, ti3name))
		error("CGATS file read error on '%s': %s",ti3name,icg->err);

	if (icg->ntables == 0 || icg->t[0].tt != tt_other || icg->t[0].oi != 0)
		error ("Input file '%s' isn't a CTI3 format file",ti3name);
	if (icg->ntables < 1)
		error ("Input file '%s' doesn't contain at least one table",ti3name);

	/* See if CIE is actually available - some sources of .TI3 don't provide it */
	if (!spec
	 && icg->find_field(icg, 0, "LAB_L") < 0
	 && icg->find_field(icg, 0, "XYZ_X") < 0) {

		if (icg->find_kword(icg, 0, "SPECTRAL_BANDS") < 0)
			error ("Neither CIE nor spectral data found in file '%s'",ti3name);

		/* Switch to using spectral information */
		if (verb)
			printf("No CIE data found, switching to spectral with standard observer & D50\n");
		spec = 1;
		illum = icxIT_D50;
		observ = icxOT_CIE_1931_2;
	}
	
	/* Figure out what sort of device it is */
	{
		int ti;

		if ((ti = icg->find_kword(icg, 0, "DEVICE_CLASS")) < 0)
			error ("Input file '%s' doesn't contain keyword DEVICE_CLASS",ti3name);

		if (strcmp(icg->t[0].kdata[ti],"DISPLAY") == 0) {
			isdisp = 1;
		}

		/* See if the CIE data has been normalised to Y = 100 */
		if ((ti = icg->find_kword(icg, 0, "NORMALIZED_TO_Y_100")) < 0
		 || strcmp(icg->t[0].kdata[ti],"NO") == 0) {
			isdnormed = 0;
		} else {
			isdnormed = 1;
		}
	}

	/* Figure out what sort of device colorspace it is */
	{
		int ti;

		if ((ti = icg->find_kword(icg, 0, "COLOR_REP")) < 0)
			error("Input file '%s' doesn't contain keyword COLOR_REPS",ti3name);

		if (strcmp(icg->t[0].kdata[ti],"CMYK_XYZ") == 0) {
			devspace = icSigCmykData;
			devchan = 4;
			isLab = 0;
			isAdditive = 0;
		} else if (strcmp(icg->t[0].kdata[ti],"CMYK_LAB") == 0) {
			devspace = icSigCmykData;
			devchan = 4;
			isLab = 1;
			isAdditive = 0;
		} else if (strcmp(icg->t[0].kdata[ti],"CMY_XYZ") == 0) {
			devspace = icSigCmyData;
			devchan = 3;
			isLab = 0;
			isAdditive = 0;
		} else if (strcmp(icg->t[0].kdata[ti],"CMY_LAB") == 0) {
			devspace = icSigCmyData;
			devchan = 3;
			isLab = 1;
			isAdditive = 0;
		} else if (strcmp(icg->t[0].kdata[ti],"RGB_XYZ") == 0) {
			devspace = icSigRgbData;
			devchan = 3;
			isLab = 0;
			isAdditive = 1;
		} else if (strcmp(icg->t[0].kdata[ti],"RGB_LAB") == 0) {
			devspace = icSigRgbData;
			devchan = 3;
			isLab = 1;
			isAdditive = 1;
		} else if (strcmp(icg->t[0].kdata[ti],"iRGB_XYZ") == 0) {
			devspace = icSigRgbData;
			devchan = 3;
			isLab = 0;
			isAdditive = 1;
		} else if (strcmp(icg->t[0].kdata[ti],"iRGB_LAB") == 0) {
			devspace = icSigRgbData;
			devchan = 3;
			isLab = 1;
			isAdditive = 1;
		/* Scanner .ti3 files: */
		} else if (strcmp(icg->t[0].kdata[ti],"XYZ_RGB") == 0) {
			devspace = icSigRgbData;
			devchan = 3;
			isLab = 0;
			isAdditive = 1;
		} else if (strcmp(icg->t[0].kdata[ti],"LAB_RGB") == 0) {
			devspace = icSigRgbData;
			devchan = 3;
			isLab = 1;
			isAdditive = 1;
#ifdef IMP_MONO
		} else if (strcmp(icg->t[0].kdata[ti],"K_XYZ") == 0) {
			devspace = icSigGrayData;
			devchan = 1;
			isLab = 0;
			isAdditive = 0;
		} else if (strcmp(icg->t[0].kdata[ti],"K_LAB") == 0) {
			devspace = icSigGrayData;
			devchan = 1;
			isLab = 1;
			isAdditive = 0;
		} else if (strcmp(icg->t[0].kdata[ti],"W_XYZ") == 0) {
			devspace = icSigGrayData;
			devchan = 1;
			isLab = 0;
			isAdditive = 1;
		} else if (strcmp(icg->t[0].kdata[ti],"W_LAB") == 0) {
			devspace = icSigGrayData;
			devchan = 1;
			isLab = 1;
			isAdditive = 1;
#endif /* IMP_MONO */

		} else 
			error("Device input file '%s' has unhandled color representation '%s'",
			                                                     ti3name, icg->t[0].kdata[ti]);
	}

	if ((npat = icg->t[0].nsets) <= 0)
		error("Input file '%s' has no sets of data",ti3name);

	if (verb) {
		fprintf(verbo,"No of test patches = %d\n",npat);
	}

	/* Allocate arrays to hold test patch input and output values */
	if ((tpat = (pval *)malloc(sizeof(pval) * npat)) == NULL)
		error("Malloc failed - tpat[]");

	/* Read in the CGATs fields */
	{
		int sidx;					/* Sample ID index */
		int sloc;					/* Sample location index (if any) */
		int ti, ci, mi, yi, ki;
		int Xi, Yi, Zi;

		if ((sidx = icg->find_field(icg, 0, "SAMPLE_ID")) < 0)
			error("Input file '%s' doesn't contain field SAMPLE_ID",ti3name);
		if (icg->t[0].ftype[sidx] != nqcs_t)
			error("Input file '%s' field SAMPLE_ID is wrong type",ti3name);

		if ((sloc = icg->find_field(icg, 0, "SAMPLE_LOC")) >= 0) {
			if (icg->t[0].ftype[sloc] != cs_t)
				error("Input file '%s' field SAMPLE_LOC is wrong type",ti3name);
		}

		if (devspace == icSigGrayData) {
			if (isAdditive) {
				if ((ci = icg->find_field(icg, 0, "GRAY_W")) < 0)
					error("Input file doesn't contain field GRAY_W");
				if (icg->t[0].ftype[ci] != r_t)
					error("Field GRAY_W is wrong type - corrupted file ?");
			} else {
				if ((ci = icg->find_field(icg, 0, "GRAY_K")) < 0)
					error("Input file doesn't contain field GRAY_K");
				if (icg->t[0].ftype[ci] != r_t)
					error("Field GRAY_K is wrong type - corrupted file ?");
			}
			mi = yi = ki = ci;

		} else if (devspace == icSigRgbData) {
			if ((ci = icg->find_field(icg, 0, "RGB_R")) < 0)
				error("Input file '%s' doesn't contain field RGB_R",ti3name);
			if (icg->t[0].ftype[ci] != r_t)
				error("Input file '%s' field RGB_R is wrong type",ti3name);
			if ((mi = icg->find_field(icg, 0, "RGB_G")) < 0)
				error("Input file '%s' doesn't contain field RGB_G",ti3name);
			if (icg->t[0].ftype[mi] != r_t)
				error("Input file '%s' field RGB_G is wrong type",ti3name);
			if ((yi = icg->find_field(icg, 0, "RGB_B")) < 0)
				error("Input file '%s' doesn't contain field RGB_B",ti3name);
			if (icg->t[0].ftype[yi] != r_t)
				error("Input file '%s' field RGB_B is wrong type",ti3name);
			ki = yi;

		} else if (devspace == icSigCmyData) {

			if ((ci = icg->find_field(icg, 0, "CMY_C")) < 0)
				error("Input file '%s' doesn't contain field CMY_C",ti3name);
			if (icg->t[0].ftype[ci] != r_t)
				error("Input file '%s' field CMY_C is wrong type",ti3name);
			if ((mi = icg->find_field(icg, 0, "CMY_M")) < 0)
				error("Input file '%s' doesn't contain field CMY_M",ti3name);
			if (icg->t[0].ftype[mi] != r_t)
				error("Input file '%s' field CMY_M is wrong type",ti3name);
			if ((yi = icg->find_field(icg, 0, "CMY_Y")) < 0)
				error("Input file '%s' doesn't contain field CMY_Y",ti3name);
			ki = yi;
		} else {	/* Assume CMYK */

			if ((ci = icg->find_field(icg, 0, "CMYK_C")) < 0)
				error("Input file '%s' doesn't contain field CMYK_C",ti3name);
			if (icg->t[0].ftype[ci] != r_t)
				error("Input file '%s' field CMYK_C is wrong type",ti3name);
			if ((mi = icg->find_field(icg, 0, "CMYK_M")) < 0)
				error("Input file '%s' doesn't contain field CMYK_M",ti3name);
			if (icg->t[0].ftype[mi] != r_t)
				error("Input file '%s' field CMYK_M is wrong type",ti3name);
			if ((yi = icg->find_field(icg, 0, "CMYK_Y")) < 0)
				error("Input file '%s' doesn't contain field CMYK_Y",ti3name);
			if (icg->t[0].ftype[yi] != r_t)
				error("Input file '%s' field CMYK_Y is wrong type",ti3name);
			if ((ki = icg->find_field(icg, 0, "CMYK_K")) < 0)
				error("Input file '%s' doesn't contain field CMYK_K",ti3name);
			if (icg->t[0].ftype[ki] != r_t)
				error("Input file '%s' field CMYK_K is wrong type",ti3name);
		}

		if (spec == 0) { 		/* Using instrument tristimulous value */

			if (isLab) {		/* Expect Lab */
				if ((Xi = icg->find_field(icg, 0, "LAB_L")) < 0)
					error("Input file '%s' doesn't contain field LAB_L",ti3name);
				if (icg->t[0].ftype[Xi] != r_t)
					error("Input file '%s' field LAB_L is wrong type",ti3name);
				if ((Yi = icg->find_field(icg, 0, "LAB_A")) < 0)
					error("Input '%s' file doesn't contain field LAB_A",ti3name);
				if (icg->t[0].ftype[Yi] != r_t)
					error("Input file '%s' field LAB_A is wrong type",ti3name);
				if ((Zi = icg->find_field(icg, 0, "LAB_B")) < 0)
					error("Input file '%s' doesn't contain field LAB_B",ti3name);
				if (icg->t[0].ftype[Zi] != r_t)
					error("Input file '%s' field LAB_B is wrong type",ti3name);

			} else { 		/* Expect XYZ */
				if ((Xi = icg->find_field(icg, 0, "XYZ_X")) < 0)
					error("Input file '%s' doesn't contain field XYZ_X",ti3name);
				if (icg->t[0].ftype[Xi] != r_t)
					error("Input file '%s' field XYZ_X is wrong type",ti3name);
				if ((Yi = icg->find_field(icg, 0, "XYZ_Y")) < 0)
					error("Input file '%s' doesn't contain field XYZ_Y",ti3name);
				if (icg->t[0].ftype[Yi] != r_t)
					error("Input file '%s' field XYZ_Y is wrong type",ti3name);
				if ((Zi = icg->find_field(icg, 0, "XYZ_Z")) < 0)
					error("Input file '%s' doesn't contain field XYZ_Z",ti3name);
				if (icg->t[0].ftype[Zi] != r_t)
					error("Input file '%s' field XYZ_Z is wrong type",ti3name);
			}

			for (i = 0; i < npat; i++) {
				strcpy(tpat[i].sid, (char *)icg->t[0].fdata[i][sidx]);
				if (sloc >= 0) 
					strcpy(tpat[i].slo, (char *)icg->t[0].fdata[i][sloc]);
				else
					strcpy(tpat[i].slo, "");
				tpat[i].p[0] = *((double *)icg->t[0].fdata[i][ci]) / 100.0;
				tpat[i].p[1] = *((double *)icg->t[0].fdata[i][mi]) / 100.0;
				tpat[i].p[2] = *((double *)icg->t[0].fdata[i][yi]) / 100.0;
				tpat[i].p[3] = *((double *)icg->t[0].fdata[i][ki]) / 100.0;
				if (tpat[i].p[0] > 1.0
				 || tpat[i].p[1] > 1.0
				 || tpat[i].p[2] > 1.0
				 || tpat[i].p[3] > 1.0) {
					error("Input file '%s' device value field value exceeds 100.0 !",ti3name);
				}
				tpat[i].v[0] = *((double *)icg->t[0].fdata[i][Xi]);
				tpat[i].v[1] = *((double *)icg->t[0].fdata[i][Yi]);
				tpat[i].v[2] = *((double *)icg->t[0].fdata[i][Zi]);
				if (!isLab && (!isdisp || isdnormed != 0)) {
					tpat[i].v[0] /= 100.0;		/* Normalise XYZ to range 0.0 - 1.0 */
					tpat[i].v[1] /= 100.0;
					tpat[i].v[2] /= 100.0;
				}
				if (!isLab) { /* Convert test patch result XYZ to PCS (D50 Lab) */
					icmXYZ2Lab(&icmD50, tpat[i].v, tpat[i].v);
				}
			}

		} else { 		/* Using spectral data */
			int ii;
			xspect sp;
			char buf[100];
			int  spi[XSPECT_MAX_BANDS];	/* CGATS indexes for each wavelength */
			xsp2cie *sp2cie;	/* Spectral conversion object */

			if ((ii = icg->find_kword(icg, 0, "SPECTRAL_BANDS")) < 0)
				error ("Input file '%s' doesn't contain keyword SPECTRAL_BANDS",ti3name);
			sp.spec_n = atoi(icg->t[0].kdata[ii]);
			if ((ii = icg->find_kword(icg, 0, "SPECTRAL_START_NM")) < 0)
				error ("Input file '%s' doesn't contain keyword SPECTRAL_START_NM",ti3name);
			sp.spec_wl_short = atof(icg->t[0].kdata[ii]);
			if ((ii = icg->find_kword(icg, 0, "SPECTRAL_END_NM")) < 0)
				error ("Input file '%s; doesn't contain keyword SPECTRAL_END_NM",ti3name);
			sp.spec_wl_long = atof(icg->t[0].kdata[ii]);
			if (!isdisp || isdnormed != 0)
				sp.norm = 100.0;
			else
				sp.norm = 1.0;

			/* Find the fields for spectral values */
			for (j = 0; j < sp.spec_n; j++) {
				int nm;
		
				/* Compute nearest integer wavelength */
				nm = (int)(sp.spec_wl_short + ((double)j/(sp.spec_n-1.0))
				            * (sp.spec_wl_long - sp.spec_wl_short) + 0.5);
				
				sprintf(buf,"SPEC_%03d",nm);

				if ((spi[j] = icg->find_field(icg, 0, buf)) < 0)
					error("Input file '%s' doesn't contain field %s",ti3name,buf);
			}

			if (isdisp) {
				illum = icxIT_none;		/* Displays are assumed to be self luminous */
			}

			/* Create a spectral conversion object */
			if ((sp2cie = new_xsp2cie(illum, illum == icxIT_none ? NULL : &cust_illum,
			                          observ, NULL, icSigLabData, icxClamp)) == NULL)
				error("Creation of spectral conversion object failed");

			if (fwacomp) {
				double nw = 0.0;		/* Number of media white patches */
				xspect mwsp;			/* Medium spectrum */
				instType itype;			/* Spectral instrument type */
				xspect insp;			/* Instrument illuminant */

				mwsp = sp;		/* Struct copy */

	 			if ((ti = icg->find_kword(icg, 0, "TARGET_INSTRUMENT")) < 0)
					error ("Input file '%s' can't find target instrument needed for FWA compensation",ti3name);

				if ((itype = inst_enum(icg->t[0].kdata[ti])) == instUnknown)
					error ("Input file '%s' unrecognised target instrument '%s'",ti3name, icg->t[0].kdata[ti]);

				if (inst_illuminant(&insp, itype) != 0)
					error ("Instrument doesn't have an FWA illuminent");

				/* Find the media white spectral reflectance */
				for (j = 0; j < mwsp.spec_n; j++)
					mwsp.spec[j] = 0.0;

				/* Compute the mean of all the media white patches */
				for (i = 0; i < npat; i++) {
					int use = 0;
	
					if (devspace == icSigGrayData) {
						if (isAdditive) {
							if (*((double *)icg->t[0].fdata[i][ci]) > (100.0 - 0.1))
								use = 1;
						} else {
							if (*((double *)icg->t[0].fdata[i][ci]) < 0.1)
								use = 1;
						}
					} else if (devspace == icSigRgbData) {
						if (*((double *)icg->t[0].fdata[i][ci]) > (100.0 - 0.1)
						 && *((double *)icg->t[0].fdata[i][mi]) > (100.0 - 0.1)
						 && *((double *)icg->t[0].fdata[i][yi]) > (100.0 - 0.1))
							use = 1;
					} else if (devspace == icSigCmyData) {
						if (*((double *)icg->t[0].fdata[i][ci]) < 0.1
						 && *((double *)icg->t[0].fdata[i][mi]) < 0.1
						 && *((double *)icg->t[0].fdata[i][yi]) < 0.1)
							use = 1;
					} else {	/* Assume CMYK */

						if (*((double *)icg->t[0].fdata[i][ci]) < 0.1
						 && *((double *)icg->t[0].fdata[i][mi]) < 0.1
						 && *((double *)icg->t[0].fdata[i][yi]) < 0.1
						 && *((double *)icg->t[0].fdata[i][ki]) < 0.1) {
							use = 1;
						}
					}

					if (use) {
						/* Read the spectral values for this patch */
						for (j = 0; j < mwsp.spec_n; j++) {
							mwsp.spec[j] += *((double *)icg->t[0].fdata[i][spi[j]]);
						}
						nw++;
					}
				}
				if (nw == 0.0) {
					warning("Input file '%s' can't find a media white patch to init FWA",ti3name);

					/* Track the maximum reflectance for any band to determine white. */
					/* This might give bogus results if there is no white patch... */
					for (i = 0; i < npat; i++) {
						for (j = 0; j < mwsp.spec_n; j++) {
							double rv = *((double *)icg->t[0].fdata[i][spi[j]]);
							if (rv > mwsp.spec[j])
								mwsp.spec[j] = rv;
						}
					}
					nw++;
				}

				for (j = 0; j < mwsp.spec_n; j++)
					mwsp.spec[j] /= nw;	/* Compute average */

				/* If we are setting a specific simulated instrument illuminant */
				if (tillum != icxIT_none) {
					tillump = &cust_tillum;
					if (tillum != icxIT_custom) {
						if (standardIlluminant(tillump, tillum, 0.0)) {
							error("simulated inst. illum. not recognised");
						}
					}
				}

				if (sp2cie->set_fwa(sp2cie, &insp, tillump, &mwsp)) 
					error ("Set FWA on sp2cie failed");

				if (verb) {
					double FWAc;
					sp2cie->get_fwa_info(sp2cie, &FWAc);
					fprintf(verbo,"FWA content = %f\n",FWAc);
				}
				
			}

			for (i = 0; i < npat; i++) {
				strcpy(tpat[i].sid, (char *)icg->t[0].fdata[i][sidx]);
				if (sloc >= 0) 
					strcpy(tpat[i].slo, (char *)icg->t[0].fdata[i][sloc]);
				else
					strcpy(tpat[i].slo, "");
				tpat[i].p[0] = *((double *)icg->t[0].fdata[i][ci]) / 100.0;
				tpat[i].p[1] = *((double *)icg->t[0].fdata[i][mi]) / 100.0;
				tpat[i].p[2] = *((double *)icg->t[0].fdata[i][yi]) / 100.0;
				tpat[i].p[3] = *((double *)icg->t[0].fdata[i][ki]) / 100.0;
				if (tpat[i].p[0] > 1.0
				 || tpat[i].p[1] > 1.0
				 || tpat[i].p[2] > 1.0
				 || tpat[i].p[3] > 1.0) {
					error("Input file '%s' device value field value exceeds 100.0 !",ti3name);
				}

				/* Read the spectral values for this patch */
				for (j = 0; j < sp.spec_n; j++) {
					sp.spec[j] = *((double *)icg->t[0].fdata[i][spi[j]]);
				}

				/* Convert it to CIE space */
				sp2cie->convert(sp2cie, tpat[i].v, &sp);
			}

			sp2cie->del(sp2cie);		/* Done with this */
		}
		/* Normalize display values to Y = 1.0 if needed */
		/* (re-norm spec derived, since observer may be different) */
		if (isdisp && (isdnormed == 0 || spec != 0)) {
			double scale = -1e6;
			double bxyz[3];

			/* Locate max Y */
			for (i = 0; i < npat; i++) {
				icmLab2XYZ(&icmD50, bxyz,  tpat[i].v);
				if (bxyz[1] > scale)
					scale = bxyz[1];
			}
			
			scale = 1.0/scale;

			/* Scale max Y to 1.0 */
			for (i = 0; i < npat; i++) {
				icmLab2XYZ(&icmD50, tpat[i].v, tpat[i].v);
				tpat[i].v[0] *= scale;
				tpat[i].v[1] *= scale;
				tpat[i].v[2] *= scale;
				icmXYZ2Lab(&icmD50, tpat[i].v, tpat[i].v);
			}
		}

		icg->del(icg);		/* Clean up */
	}	/* End of reading in CGATs file */

	/* - - - - - - - - - - */
	/* Check the forward profile accuracy against the data points */
	{
		double merr = 0.0;		/* Max */
		double aerr = 0.0;		/* Avg */
		double rerr = 0.0;		/* RMS */
		double nsamps = 0.0;
		int inn, outn;			/* Chanells for input and output spaces */

		if (dovrml) {
			wrl = start_vrml(out_name, doaxes);
			start_line_set(wrl);
		}

		/* Open up the file for reading */
		if ((rd_fp = new_icmFileStd_name(iccname,"r")) == NULL)
			error("Write: Can't open file '%s'",iccname);

		if ((rd_icco = new_icc()) == NULL)
			error("Read: Creation of ICC object failed");

		/* Read the header and tag list */
		if ((rv = rd_icco->read(rd_icco,rd_fp,0)) != 0)
			error("Read: %d, %s",rv,rd_icco->err);

		/* Get the Fwd table, absolute with Lab override */
		if ((luo = rd_icco->get_luobj(rd_icco, icmFwd, intent,
		                              icSigLabData, icmLuOrdNorm)) == NULL) {
			error("%d, %s",rd_icco->errc, rd_icco->err);
		}

		/* Get details of conversion (Arguments may be NULL if info not needed) */
		luo->spaces(luo, NULL, &inn, NULL, &outn, NULL, NULL, NULL, NULL, NULL);

		for (i = 0; i < npat; i++) {
			double out[3];
			double mxd;

			/* Lookup the patch value in the profile */
			if (luo->lookup(luo, out, tpat[i].p) > 1)
				error("%d, %s",rd_icco->errc,rd_icco->err);

			if (verb > 1) {
				printf("[%f] %s%s%s: %s -> %f %f %f should be %f %f %f\n",
				       cie2k ? icmCIE2K(tpat[i].v, out) : 
				               cie94 ? icmCIE94(tpat[i].v, out) : icmLabDE(tpat[i].v, out),
				       tpat[i].sid,
				       tpat[i].slo[0] != '\000' ? " @ " : "",
				       tpat[i].slo,
				       icmPdv(devchan, tpat[i].p),
				       out[0],out[1],out[2],
				       tpat[i].v[0],tpat[i].v[1],tpat[i].v[2]);
			}
			if (dovrml) {
				if (dominl && icmLabDE(tpat[i].v, out) < 0.5) {
					double cent[3], vec[3], vlen;
					double p1[3], p2[3];

					/* Compute center */
					icmAdd3(cent, tpat[i].v, out);
					icmScale3(cent, cent, 0.5);
					if ((vlen = icmLabDE(tpat[i].v, out)) < 1e-6) {
						vec[0] = 0.25; vec[1] = 0.0; vec[2] = 0.0;
					} else {
						icmSub3(vec, tpat[i].v, out);
						icmScale3(vec, vec, 0.25/vlen);
					}
					icmSub3(p1, cent, vec);
					icmAdd3(p2, cent, vec);
					add_vertex(wrl, p1);
					add_vertex(wrl, p2);
				} else {
					add_vertex(wrl, tpat[i].v);
					add_vertex(wrl, out);
				}
			}

			/* Check the result */
			if (cie2k)
				mxd = icmCIE2K(tpat[i].v, out);
			else if (cie94)
				mxd = icmCIE94(tpat[i].v, out);
			else
				mxd = icmLabDE(tpat[i].v, out);

			aerr += mxd;
			rerr += mxd * mxd;

			nsamps++;
			if (mxd > merr)
				merr = mxd;

		}
		if (dovrml) {
			if (dodecol)
				make_de_lines(wrl);
			else
				make_lines(wrl, 2);
			end_vrml(wrl);
		}
		printf("Profile check complete, errors%s: max. = %f, avg. = %f, RMS = %f\n",
            cie2k ? "(CIEDE2000)" : cie94 ? " (CIE94)" : "", merr, aerr/nsamps, sqrt(rerr/nsamps));

		/* ------------------------------- */
		/* If we want sort by target value */
		if (ddevv) {
			double cieval[3];

			/* Lookup the CIE value of the target */
			if (luo->lookup(luo, cieval, devval) > 1)
				error("%d, %s",rd_icco->errc,rd_icco->err);

			/* Compute deltas to target value. */
			for (i = 0; i < npat; i++) {
				if (cie2k)
					tpat[i].dv = icmCIE2K(tpat[i].v, cieval);
				else if (cie94)
					tpat[i].dv = icmCIE94(tpat[i].v, cieval);
				else
					tpat[i].dv = icmLabDE(tpat[i].v, cieval);

				tpat[i].dp = 0.0;
				for (j = 0; j < inn; j++) {
					double tt;
					tt = tpat[i].p[j] - devval[j];
					tpat[i].dp += tt * tt;
				}
				tpat[i].dp = sqrt(tpat[i].dp);
			}

			if (sortbypcs) {
				/* Sort by pcs delta */
#define HEAP_COMPARE(A,B) (A.dv < B.dv)
				HEAPSORT(pval, tpat, npat);
#undef HEAP_COMPARE
			} else {
				/* Sort by device delta */
#define HEAP_COMPARE(A,B) (A.dp < B.dp)
				HEAPSORT(pval, tpat, npat);
#undef HEAP_COMPARE
			}

			printf("Target point:\n");
			if (devspace == icSigCmykData) {
				printf("%f %f %f %f -> %f %f %f\n",
				       devval[0],devval[1],devval[2],devval[3],
				       cieval[0],cieval[1],cieval[2]);
			} else {	/* Assume RGB/CMY */
				printf("%f %f %f -> %f %f %f\n",
				       devval[0],devval[1],devval[2],
				       cieval[0],cieval[1],cieval[2]);
			}
			printf("\n");
	
			for (i = 0; i < npat; i++) {
				if (devspace == icSigCmykData) {
					printf("%s: %f %f %f %f [%f] -> %f %f %f [%f]\n",
					       tpat[i].sid,
					       tpat[i].p[0],tpat[i].p[1],tpat[i].p[2],tpat[i].p[3],
					       tpat[i].dp,
					       tpat[i].v[0],tpat[i].v[1],tpat[i].v[2],
					       tpat[i].dv);
				} else {	/* Assume RGB/CMY */
					printf("%s: %f %f %f [%f] -> %f %f %f [%f]\n",
					       tpat[i].sid,
					       tpat[i].p[0],tpat[i].p[1],tpat[i].p[2],
					       tpat[i].dp,
					       tpat[i].v[0],tpat[i].v[1],tpat[i].v[2],
					       tpat[i].dv);
				}
			}
		}

		/* Done with lookup object */
		luo->del(luo);

		/* Close the file */
		rd_icco->del(rd_icco);
		rd_fp->del(rd_fp);
	}

	return 0;
}


/* ------------------------------------------------ */
/* Some simple functions to do basix VRML work */
/* !!! Should change to plot/vrml lib !!! */

#define GAMUT_LCENT 50.0
static int npoints = 0;
static int paloc = 0;
static struct { double pp[3]; } *pary;

static void Lab2RGB(double *out, double *in);
static void DE2RGB(double *out, double in);

FILE *start_vrml(char *name, int doaxes) {
	FILE *wrl;

	/* Define the axis boxes */
	struct {
		double x, y, z;			/* Box center */
		double wx, wy, wz;		/* Box size */
		double r, g, b;			/* Box color */
	} axes[5] = {
		{ 0, 0,   50-GAMUT_LCENT, 2, 2, 100, .7, .7, .7 },	/* L axis */
		{ 50, 0,  0-GAMUT_LCENT,  100, 2, 2,  1,  0,  0 },	/* +a (red) axis */
		{ 0, -50, 0-GAMUT_LCENT,  2, 100, 2,  0,  0,  1 },	/* -b (blue) axis */
		{ -50, 0, 0-GAMUT_LCENT,  100, 2, 2,  0,  1,  0 },	/* -a (green) axis */
		{ 0,  50, 0-GAMUT_LCENT,  2, 100, 2,  1,  1,  0 },	/* +b (yellow) axis */
	};

	/* Define the labels */
	struct {
		double x, y, z;
		double size;
		char *string;
		double r, g, b;
	} labels[6] = {
		{ -2, 2, -GAMUT_LCENT + 100 + 10, 10, "+L*",  .7, .7, .7 },	/* Top of L axis */
		{ -2, 2, -GAMUT_LCENT - 10,      10, "0",    .7, .7, .7 },	/* Bottom of L axis */
		{ 100 + 5, -3,  0-GAMUT_LCENT,  10, "+a*",  1,  0,  0 },	/* +a (red) axis */
		{ -5, -100 - 10, 0-GAMUT_LCENT,  10, "-b*",  0,  0,  1 },	/* -b (blue) axis */
		{ -100 - 15, -3, 0-GAMUT_LCENT,  10, "-a*",  0,  0,  1 },	/* -a (green) axis */
		{ -5,  100 + 5, 0-GAMUT_LCENT,  10, "+b*",  1,  1,  0 },	/* +b (yellow) axis */
	};

	if ((wrl = fopen(name,"w")) == NULL)
		error("Error opening VRML file '%s'\n",name);

	npoints = 0;

	fprintf(wrl,"#VRML V2.0 utf8\n");
	fprintf(wrl,"\n");
	fprintf(wrl,"# Created by the Argyll CMS\n");
	fprintf(wrl,"Transform {\n");
	fprintf(wrl,"children [\n");
	fprintf(wrl,"	NavigationInfo {\n");
	fprintf(wrl,"		type \"EXAMINE\"        # It's an object we examine\n");
	fprintf(wrl,"	} # We'll add our own light\n");
	fprintf(wrl,"\n");
	fprintf(wrl,"    DirectionalLight {\n");
	fprintf(wrl,"        direction 0 0 -1      # Light illuminating the scene\n");
	fprintf(wrl,"        direction 0 -1 0      # Light illuminating the scene\n");
	fprintf(wrl,"    }\n");
	fprintf(wrl,"\n");
	fprintf(wrl,"    Viewpoint {\n");
	fprintf(wrl,"        position 0 0 340      # Position we view from\n");
	fprintf(wrl,"    }\n");
	fprintf(wrl,"\n");
	if (doaxes != 0) {
		int n;
		fprintf(wrl,"    # Lab axes as boxes:\n");
		for (n = 0; n < 5; n++) {
			fprintf(wrl,"    Transform { translation %f %f %f\n", axes[n].x, axes[n].y, axes[n].z);
			fprintf(wrl,"      children [\n");
			fprintf(wrl,"        Shape{\n");
			fprintf(wrl,"          geometry Box { size %f %f %f }\n",
			                       axes[n].wx, axes[n].wy, axes[n].wz);
			fprintf(wrl,"          appearance Appearance { material Material ");
			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[n].r, axes[n].g, axes[n].b);
			fprintf(wrl,"        }\n");
			fprintf(wrl,"      ]\n");
			fprintf(wrl,"    }\n");
		}
		fprintf(wrl,"    # Axes identification:\n");
		for (n = 0; n < 6; n++) {
			fprintf(wrl,"    Transform { translation %f %f %f\n", labels[n].x, labels[n].y, labels[n].z);
			fprintf(wrl,"      children [\n");
			fprintf(wrl,"        Shape{\n");
			fprintf(wrl,"          geometry Text { string [\"%s\"]\n",labels[n].string);
			fprintf(wrl,"            fontStyle FontStyle { family \"SANS\" style \"BOLD\" size %f }\n",
			                                  labels[n].size);
			fprintf(wrl,"                        }\n");
			fprintf(wrl,"          appearance Appearance { material Material ");
			fprintf(wrl,"{ diffuseColor %f %f %f} }\n", labels[n].r, labels[n].g, labels[n].b);
			fprintf(wrl,"        }\n");
			fprintf(wrl,"      ]\n");
			fprintf(wrl,"    }\n");
		}
		fprintf(wrl,"\n");
	}

	return wrl;
}

void
start_line_set(FILE *wrl) {

	fprintf(wrl,"\n");
	fprintf(wrl,"Shape {\n");
	fprintf(wrl,"  geometry IndexedLineSet { \n");
	fprintf(wrl,"    coord Coordinate { \n");
	fprintf(wrl,"	   point [\n");
}

void add_vertex(FILE *wrl, double pp[3]) {

	fprintf(wrl,"%f %f %f,\n",pp[1], pp[2], pp[0]-GAMUT_LCENT);
	
	if (paloc < (npoints+1)) {
		paloc = (paloc + 10) * 2;
		if (pary == NULL)
			pary = malloc(paloc * 3 * sizeof(double));
		else
			pary = realloc(pary, paloc * 3 * sizeof(double));

		if (pary == NULL)
			error ("Malloc failed");
	}
	pary[npoints].pp[0] = pp[0];
	pary[npoints].pp[1] = pp[1];
	pary[npoints].pp[2] = pp[2];
	npoints++;
}


void make_lines(FILE *wrl, int ppset) {
	int i, j;

	fprintf(wrl,"      ]\n");
	fprintf(wrl,"    }\n");
	fprintf(wrl,"  coordIndex [\n");

	for (i = 0; i < npoints;) {
		for (j = 0; j < ppset; j++, i++) {
			fprintf(wrl,"%d, ", i);
		}
		fprintf(wrl,"-1,\n");
	}
	fprintf(wrl,"    ]\n");

	/* Color */
	fprintf(wrl,"            colorPerVertex TRUE\n");
	fprintf(wrl,"            color Color {\n");
	fprintf(wrl,"              color [			# RGB colors of each vertex\n");

	for (i = 0; i < npoints; i++) {
		double rgb[3], Lab[3];
		Lab[0] = pary[i].pp[0];
		Lab[1] = pary[i].pp[1];
		Lab[2] = pary[i].pp[2];
		Lab2RGB(rgb, Lab);
		fprintf(wrl,"                %f %f %f,\n", rgb[0], rgb[1], rgb[2]);
	}
	fprintf(wrl,"              ] \n");
	fprintf(wrl,"            }\n");
	/* End color */

	fprintf(wrl,"  }\n");
	fprintf(wrl,"} # end shape\n");
}

/* Assume 2 ppset, and make line color prop to length */
void make_de_lines(FILE *wrl) {
	int i, j;

	fprintf(wrl,"      ]\n");
	fprintf(wrl,"    }\n");
	fprintf(wrl,"  coordIndex [\n");

	for (i = 0; i < npoints;) {
		for (j = 0; j < 2; j++, i++) {
			fprintf(wrl,"%d, ", i);
		}
		fprintf(wrl,"-1,\n");
	}
	fprintf(wrl,"    ]\n");

	/* Color */
	fprintf(wrl,"            colorPerVertex TRUE\n");
	fprintf(wrl,"            color Color {\n");
	fprintf(wrl,"              color [			# RGB colors of each vertex\n");

	for (i = 0; i < npoints; i++) {
		double rgb[3], ss;
		for (ss = 0.0, j = 0; j < 3; j++) {
			double tt = (pary[i & ~1].pp[j] - pary[i | 1].pp[j]);
			ss += tt * tt;
		}
		ss = sqrt(ss);
		DE2RGB(rgb, ss);
		fprintf(wrl,"                %f %f %f,\n", rgb[0], rgb[1], rgb[2]);
	}
	fprintf(wrl,"              ] \n");
	fprintf(wrl,"            }\n");
	/* End color */

	fprintf(wrl,"  }\n");
	fprintf(wrl,"} # end shape\n");
}

void end_vrml(FILE *wrl) {

	fprintf(wrl,"\n");
	fprintf(wrl,"  ] # end of children for world\n");
	fprintf(wrl,"}\n");

	if (fclose(wrl) != 0)
		error("Error closing VRML file\n");
}


/* Convert a gamut Lab value to an RGB value for display purposes */
static void
Lab2RGB(double *out, double *in) {
	double L = in[0], a = in[1], b = in[2];
	double x,y,z,fx,fy,fz;
	double R, G, B;

	/* Scale so that black is visible */
	L = L * (100 - 40.0)/100.0 + 40.0;

	/* First convert to XYZ using D50 white point */
	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;

	x *= 0.9642;	/* Multiply by white point, D50 */
	y *= 1.0;
	z *= 0.8249;

	/* Now convert to sRGB values */
	R = x * 3.2410  + y * -1.5374 + z * -0.4986;
	G = x * -0.9692 + y * 1.8760  + z * 0.0416;
	B = x * 0.0556  + y * -0.2040 + z * 1.0570;

	if (R < 0.0)
		R = 0.0;
	else if (R > 1.0)
		R = 1.0;

	if (G < 0.0)
		G = 0.0;
	else if (G > 1.0)
		G = 1.0;

	if (B < 0.0)
		B = 0.0;
	else if (B > 1.0)
		B = 1.0;

	R = pow(R, 1.0/2.2);
	G = pow(G, 1.0/2.2);
	B = pow(B, 1.0/2.2);

	out[0] = R;
	out[1] = G;
	out[2] = B;
}

/* Convert a delta E value into a signal color: */
static void
DE2RGB(double *out, double in) {
	struct {
		double de;
		double r, g, b;
	} range[6] = {
		{ 10.0, 1, 1, 0 },		/* yellow */
		{ 4.0,  1, 0, 0 },		/* red */
		{ 2.0, 1, 0, 1 },		/* magenta */
		{ 1.0, 0, 0, 1 },		/* blue */
		{ 0.5, 0, 1, 1 },		/* cyan */
		{ 0.0, 0, 1, 0 }		/* green */
	};
	int i;
	double bl;

	/* Locate the range we're in */
	if (in > range[0].de) {
		out[0] = range[0].r;
		out[1] = range[0].g;
		out[2] = range[0].b;
	} else {
		for (i = 0; i < 5; i++) {
			if (in <= range[i].de && in >= range[i+1].de)
				break;
		}
		bl = (in - range[i+1].de)/(range[i].de - range[i+1].de);
		out[0] = bl * range[i].r + (1.0 - bl) * range[i+1].r;
		out[1] = bl * range[i].g + (1.0 - bl) * range[i+1].g;
		out[2] = bl * range[i].b + (1.0 - bl) * range[i+1].b;
	}
}