1 
2 /*
3  * Argyll Color Correction System
4  * Color Device profile checker.
5  *
6  * Author: Graeme W. Gill
7  * Date:   15/7/2001
8  *
9  * Copyright 2001 - 2005 Graeme W. Gill
10  * All rights reserved.
11  *
12  * This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
13  * see the License.txt file for licencing details.
14  */
15 
16 /*
17  * This program takes in the .ti3 scattered test chart
18  * points, and checks them against an ICC profile.
19  * forward ICC device profile.
20  */
21 
22 /*
23  * TTBD:
24  *		Switch to generic colorant read code rather than Grey/RGB/CMYK,
25  *		and allow checking ICC profiles > 4 colors
26  */
27 
28 #undef DEBUG
29 
30 #undef HACK					/* Print per patch XYZ differences */
31 
32 #define IMP_MONO			/* Turn on development code */
33 
34 #define verbo stdout
35 
36 #include <stdio.h>
37 #include <string.h>
38 #include <math.h>
39 #include <ctype.h>
40 #include "copyright.h"
41 #include "aconfig.h"
42 #include "numlib.h"
43 #include "cgats.h"
44 #include "xicc.h"
45 #include "insttypes.h"
46 #include "sort.h"
47 #include "plot.h"
48 #include "vrml.h"
49 #include "ui.h"
50 
51 void
usage(void)52 usage(void) {
53 	fprintf(stderr,"Check accuracy of ICC profile, Version %s\n",ARGYLL_VERSION_STR);
54 	fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
55 	fprintf(stderr,"usage: profcheck [-options] data.ti3 iccprofile.icm\n");
56 	fprintf(stderr," -v [level]      Verbosity level (default 1), 2 to print each DE\n");
57 	fprintf(stderr," -c              Show CIE94 delta E values\n");
58 	fprintf(stderr," -k              Show CIEDE2000 delta E values\n");
59 	fprintf(stderr," -w              create %s visualisation (iccprofile%s)\n",vrml_format(),vrml_ext());
60 	fprintf(stderr," -x              Use %s axes\n",vrml_format());
61 	fprintf(stderr," -m              Make %s lines a minimum of 0.5\n",vrml_format());
62 	fprintf(stderr," -e              Color vectors acording to delta E\n");
63 	fprintf(stderr," -h              Plot a histogram of delta E's\n");
64 	fprintf(stderr," -s              Sort output by delta E\n");
65 	fprintf(stderr," -P N.NN         Create a pruned .ti3 with points less or equal to N.NN delta E\n");
66 	fprintf(stderr," -d devval1,deval2,devvalN\n");
67 	fprintf(stderr,"                 Specify a device value to sort against\n");
68 	fprintf(stderr," -p              Sort device value by PCS (Lab) target\n");
69 	fprintf(stderr," -f [illum]      Use Fluorescent Whitening Agent compensation [opt. simulated inst. illum.:\n");
70 	fprintf(stderr,"                  M0, M1, M2, A, C, D50 (def.), D50M2, D65, F5, F8, F10 or file.sp]\n");
71 	fprintf(stderr," -i illum        Choose illuminant for computation of CIE XYZ from spectral data & FWA:\n");
72 	fprintf(stderr,"                  A, C, D50 (def.), D50M2, D65, F5, F8, F10 or file.sp\n");
73 	fprintf(stderr," -o observ       Choose CIE Observer for spectral data:\n");
74 	fprintf(stderr,"                  1931_2 (def), 1964_10, S&B 1955_2, shaw, J&V 1978_2\n");
75 	fprintf(stderr," -I intent       r = relative colorimetric, a = absolute (default)\n");
76 	fprintf(stderr," data.ti3        Test data file\n");
77 	fprintf(stderr," iccprofile.icm  Profile to check against\n");
78 	exit(1);
79 	}
80 
81 static void DE2RGB(double *out, double in);
82 
83 /* Patch value type */
84 typedef struct {
85 	char sid[50];		/* sample id */
86 	char slo[50];		/* sample location, "" if not known */
87 	double p[MAX_CHAN];	/* Device value */
88 	double v[3];		/* CIE value */
89 
90 	double pv[3];		/* Profile CIE value */
91 	double de;			/* Delta E to profile CIE value */
92 
93 	double dp;			/* Delta p[] from target device value */
94 	double dv;			/* Delta E from CIE value */
95 } pval;
96 
97 /* Histogram bin type */
98 typedef struct {
99 	int count;			/* Raw count */
100 	double val;			/* Normalized value */
101 	double min, max;	/* Bin range */
102 } hbin;
103 
main(int argc,char * argv[])104 int main(int argc, char *argv[])
105 {
106 	int fa,nfa;				/* current argument we're looking at */
107 	int verb = 0;
108 	int cie94 = 0;
109 	int cie2k = 0;
110 	int dovrml = 0;
111 	int dominl = 0;
112 	int doaxes = 0;
113 	int dodecol = 0;
114 	char ti3name[MAXNAMEL+1] = { 0 };	/* Input cgats file base name */
115 	cgats *icg;				/* input cgats structure */
116 	char iccname[MAXNAMEL+1] = { 0 };	/* Input icc file base name */
117 	icmFile *rd_fp;
118 	icRenderingIntent intent = icAbsoluteColorimetric;
119 	icc *rd_icco;
120 	icmLuBase *luo;
121 	char out_name[MAXNAMEL+1], *xl;		/* VRML/X3D name */
122 	vrml *wrl = NULL;
123 
124 	int fwacomp = 0;			/* FWA compensation */
125 	int isemis = 0;				/* nz if this has emissive reference data */
126 	int isdisp = 0;				/* nz if this is a display device, 0 if output */
127 	int isdnormed = 0;      	/* Has display data been normalised to 100 ? */
128 	int spec = 0;				/* Use spectral data flag */
129 	icxIllumeType tillum = icxIT_none;	/* Target/simulated instrument illuminant */
130 	xspect cust_tillum, *tillump = NULL; /* Custom target/simulated illumination spectrum */
131 	icxIllumeType illum = icxIT_none;	/* Spectral defaults */
132 	xspect cust_illum;			/* Custom illumination spectrum */
133 	icxObserverType observ = icxOT_none;
134 
135 	int sortbyde = 0;			/* Sort by delta E */
136 
137 	int ddevv = 0;				/* Do device value sort */
138 	double devval[MAX_CHAN];	/* device value to sort on */
139 	int sortbypcs = 0;			/* Sort by PCS error to device target */
140 
141 	int dohisto = 0;			/* Plot histogram of delta E's */
142 	double prune = 0.0;			/* If > 0.0, created a pruned .ti3 file */
143 
144 	int npat;					/* Number of patches */
145 	pval *tpat;					/* Patch input values */
146 	pval **stpat;				/* Pointers to internal sorted tpat[] */
147 	int i, j, rv = 0;
148 	icColorSpaceSignature devspace = 0;	/* The device colorspace */
149 	int isAdditive = 0;			/* 0 if subtractive, 1 if additive colorspace */
150 	int isLab = 0;				/* 0 if input is XYZ, 1 if input is Lab */
151 	int devchan = 0;			/* Number of device chanels */
152 
153 #if defined(__IBMC__)
154 	_control87(EM_UNDERFLOW, EM_UNDERFLOW);
155 	_control87(EM_OVERFLOW, EM_OVERFLOW);
156 #endif
157 
158 	error_program = "profcheck";
159 
160 	if (argc <= 1)
161 		usage();
162 
163 	/* Process the arguments */
164 	for(fa = 1;fa < argc;fa++) {
165 
166 		nfa = fa;					/* skip to nfa if next argument is used */
167 		if (argv[fa][0] == '-') {	/* Look for any flags */
168 			char *na = NULL;		/* next argument after flag, null if none */
169 
170 			if (argv[fa][2] != '\000')
171 				na = &argv[fa][2];		/* next is directly after flag */
172 			else {
173 				if ((fa+1) < argc) {
174 					if (argv[fa+1][0] != '-') {
175 						nfa = fa + 1;
176 						na = argv[nfa];		/* next is seperate non-flag argument */
177 					}
178 				}
179 			}
180 
181 			if (argv[fa][1] == '?')
182 				usage();
183 
184 			/* Verbosity */
185 			else if (argv[fa][1] == 'v') {
186 				verb = 1;
187 				if (na != NULL && isdigit(na[0])) {
188 					verb = atoi(na);
189 				}
190 			}
191 
192 			/* VRML/X3D */
193 			else if (argv[fa][1] == 'w')
194 				dovrml = 1;
195 
196 			/* Minimum line length */
197 			else if (argv[fa][1] == 'm')
198 				dominl = 1;
199 
200 			/* Axes */
201 			else if (argv[fa][1] == 'x')
202 				doaxes = 1;
203 
204 			/* Delta E coloring */
205 			else if (argv[fa][1] == 'e')
206 				dodecol = 1;
207 
208 			else if (argv[fa][1] == 'c') {
209 				cie94 = 1;
210 				cie2k = 0;
211 			}
212 
213 			else if (argv[fa][1] == 'k') {
214 				cie94 = 0;
215 				cie2k = 1;
216 			}
217 
218 			/* Plot histogram */
219 			else if (argv[fa][1] == 'h')
220 				dohisto = 1;
221 
222 			/* Sort by delta E */
223 			else if (argv[fa][1] == 's')
224 				sortbyde = 1;
225 
226 			/* Create a pruned .ti3 */
227 			else if (argv[fa][1] == 'P') {
228 				if (na == NULL) usage();
229 				fa = nfa;
230 				prune = atof(na);
231 				if (prune <= 0.0)
232 					usage();
233 			}
234 
235 			/* Device sort value */
236 			else if (argv[fa][1] == 'd') {
237 				char *tp, buf[200];
238 				int ndv;
239 				if (na == NULL) usage();
240 				fa = nfa;
241 
242 				ddevv = 1;
243 				strcpy(buf, na);
244 
245 				/* Replace ',' with '\000' */
246 				for (ndv = 1,tp = buf; *tp != '\000'; tp++) {
247 					if (*tp == ',') {
248 						*tp = '\000';
249 						ndv++;
250 					}
251 				}
252 				if (ndv >= MAX_CHAN)
253 					ndv = MAX_CHAN;
254 
255 				for (tp = buf, i = 0; i < ndv; i++, tp += strlen(tp) + 1) {
256 					devval[i] = atof(tp);
257 				}
258 			}
259 
260 			else if (argv[fa][1] == 'p')
261 				sortbypcs = 1;
262 
263 			/* FWA compensation */
264 			else if (argv[fa][1] == 'f') {
265 				fwacomp = 1;
266 
267 				if (na != NULL) {	/* Argument is present - target/simulated instr. illum. */
268 					fa = nfa;
269 					if (strcmp(na, "A") == 0
270 					 || strcmp(na, "M0") == 0) {
271 						spec = 1;
272 						tillum = icxIT_A;
273 					} else if (strcmp(na, "C") == 0) {
274 						spec = 1;
275 						tillum = icxIT_C;
276 					} else if (strcmp(na, "D50") == 0
277 					        || strcmp(na, "M1") == 0) {
278 						spec = 1;
279 						tillum = icxIT_D50;
280 					} else if (strcmp(na, "D50M2") == 0
281 					        || strcmp(na, "M2") == 0) {
282 						spec = 1;
283 						tillum = icxIT_D50M2;
284 					} else if (strcmp(na, "D65") == 0) {
285 						spec = 1;
286 						tillum = icxIT_D65;
287 					} else if (strcmp(na, "F5") == 0) {
288 						spec = 1;
289 						tillum = icxIT_F5;
290 					} else if (strcmp(na, "F8") == 0) {
291 						spec = 1;
292 						tillum = icxIT_F8;
293 					} else if (strcmp(na, "F10") == 0) {
294 						spec = 1;
295 						tillum = icxIT_F10;
296 					} else {	/* Assume it's a filename */
297 						inst_meas_type mt;
298 
299 						spec = 1;
300 						tillum = icxIT_custom;
301 						if (read_xspect(&cust_tillum, &mt, na) != 0)
302 							usage();
303 
304 						if (mt != inst_mrt_none
305 						 && mt != inst_mrt_emission
306 						 && mt != inst_mrt_ambient
307 						 && mt != inst_mrt_emission_flash
308 						 && mt != inst_mrt_ambient_flash) {
309 							error("Target illuminant '%s' is wrong measurement type",na);
310 						}
311 					}
312 				}
313 			}
314 			/* Spectral Illuminant type */
315 			else if (argv[fa][1] == 'i') {
316 				if (na == NULL) usage();
317 				fa = nfa;
318 				if (strcmp(na, "A") == 0) {
319 					spec = 1;
320 					illum = icxIT_A;
321 				} else if (strcmp(na, "C") == 0) {
322 					spec = 1;
323 					illum = icxIT_C;
324 				} else if (strcmp(na, "D50") == 0) {
325 					spec = 1;
326 					illum = icxIT_D50;
327 				} else if (strcmp(na, "D50M2") == 0) {
328 					spec = 1;
329 					illum = icxIT_D50M2;
330 				} else if (strcmp(na, "D65") == 0) {
331 					spec = 1;
332 					illum = icxIT_D65;
333 				} else if (strcmp(na, "F5") == 0) {
334 					spec = 1;
335 					illum = icxIT_F5;
336 				} else if (strcmp(na, "F8") == 0) {
337 					spec = 1;
338 					illum = icxIT_F8;
339 				} else if (strcmp(na, "F10") == 0) {
340 					spec = 1;
341 					illum = icxIT_F10;
342 				} else {	/* Assume it's a filename */
343 					inst_meas_type mt;
344 
345 					spec = 1;
346 					illum = icxIT_custom;
347 					if (read_xspect(&cust_illum, &mt, na) != 0)
348 						usage();
349 
350 					if (mt != inst_mrt_none
351 					 && mt != inst_mrt_emission
352 					 && mt != inst_mrt_ambient
353 					 && mt != inst_mrt_emission_flash
354 					 && mt != inst_mrt_ambient_flash) {
355 						error("CIE illuminant '%s' is wrong measurement type",na);
356 					}
357 				}
358 			}
359 
360 			/* Spectral Observer type */
361 			else if (argv[fa][1] == 'o') {
362 				if (na == NULL) usage();
363 				fa = nfa;
364 				if (strcmp(na, "1931_2") == 0) {			/* Classic 2 degree */
365 					spec = 1;
366 					observ = icxOT_CIE_1931_2;
367 				} else if (strcmp(na, "1964_10") == 0) {	/* Classic 10 degree */
368 					spec = 1;
369 					observ = icxOT_CIE_1964_10;
370 				} else if (strcmp(na, "1955_2") == 0) {		/* Stiles and Burch 1955 2 degree */
371 					spec = 1;
372 					observ = icxOT_Stiles_Burch_2;
373 				} else if (strcmp(na, "1978_2") == 0) {		/* Judd and Voss 1978 2 degree */
374 					spec = 1;
375 					observ = icxOT_Judd_Voss_2;
376 				} else if (strcmp(na, "shaw") == 0) {		/* Shaw and Fairchilds 1997 2 degree */
377 					spec = 1;
378 					observ = icxOT_Shaw_Fairchild_2;
379 				} else
380 					usage();
381 			}
382 
383 			/* Intent (only applies to ICC profile) */
384 			else if (argv[fa][1] == 'I') {
385 				if (na == NULL) usage();
386 				fa = nfa;
387     			switch (na[0]) {
388 					case 'r':
389 						intent = icRelativeColorimetric;
390 						break;
391 					case 'a':
392 						intent = icAbsoluteColorimetric;
393 						break;
394 					default:
395 						usage();
396 				}
397 			}
398 
399 			else
400 				usage();
401 		} else
402 			break;
403 	}
404 
405 	/* Get the file name arguments */
406 	if (fa >= argc || argv[fa][0] == '-') usage();
407 	strncpy(ti3name,argv[fa++],MAXNAMEL); ti3name[MAXNAMEL] = '\000';
408 
409 	if (fa >= argc || argv[fa][0] == '-') usage();
410 	strncpy(iccname,argv[fa++],MAXNAMEL); iccname[MAXNAMEL] = '\000';
411 
412 	strncpy(out_name,iccname,MAXNAMEL-4); out_name[MAXNAMEL-4] = '\000';
413 	if ((xl = strrchr(out_name, '.')) == NULL)	/* Figure where extention is */
414 		xl = out_name + strlen(out_name);
415 	xl[0] = '\000';			/* Remove extension */
416 
417 	if (fwacomp && spec == 0)
418 		error ("FWA compensation only works when viewer and/or illuminant selected");
419 
420 	/* Open and look at the .ti3 profile patches file */
421 	icg = new_cgats();			/* Create a CGATS structure */
422 	icg->add_other(icg, "CTI3"); 	/* our special input type is Calibration Target Information 3 */
423 
424 	if (icg->read_name(icg, ti3name))
425 		error("CGATS file read error on '%s': %s",ti3name,icg->err);
426 
427 	if (icg->ntables == 0 || icg->t[0].tt != tt_other || icg->t[0].oi != 0)
428 		error ("Input file '%s' isn't a CTI3 format file",ti3name);
429 	if (icg->ntables < 1)
430 		error ("Input file '%s' doesn't contain at least one table",ti3name);
431 
432 	/* Figure out what sort of device it is */
433 	{
434 		int ti;
435 
436 		if ((ti = icg->find_kword(icg, 0, "DEVICE_CLASS")) < 0)
437 			error ("Input file '%s' doesn't contain keyword DEVICE_CLASS",ti3name);
438 
439 		if (strcmp(icg->t[0].kdata[ti],"EMISINPUT") == 0) {
440 			isemis = 1;
441 
442 		} else if (strcmp(icg->t[0].kdata[ti],"DISPLAY") == 0) {
443 			isemis = 1;
444 			isdisp = 1;
445 		}
446 
447 		/* See if the CIE data has been normalised to Y = 100 */
448 		if ((ti = icg->find_kword(icg, 0, "NORMALIZED_TO_Y_100")) < 0
449 		 || strcmp(icg->t[0].kdata[ti],"NO") == 0) {
450 			isdnormed = 0;
451 		} else {
452 			isdnormed = 1;
453 		}
454 	}
455 
456 	if (isemis && illum != icxIT_none)
457 		warning("-i illuminant ignored for emissive reference type");
458 
459 	if (isemis & fwacomp) {
460 			warning("-f FWA compensation ignored for emissive reference type");
461 		fwacomp = 0;
462 		tillum = icxIT_none;
463 	}
464 
465 	/* Set defaults */
466 	if (illum == icxIT_none)
467 		illum = icxIT_D50;
468 
469 	if (observ == icxOT_none)
470 		observ = icxOT_CIE_1931_2;
471 
472 	/* See if CIE is actually available - some sources of .TI3 don't provide it */
473 	if (!spec
474 	 && icg->find_field(icg, 0, "LAB_L") < 0
475 	 && icg->find_field(icg, 0, "XYZ_X") < 0) {
476 
477 		if (icg->find_kword(icg, 0, "SPECTRAL_BANDS") < 0)
478 			error ("Neither CIE nor spectral data found in file '%s'",ti3name);
479 
480 		/* Switch to using spectral information */
481 		if (verb)
482 			printf("No CIE data found, switching to spectral with standard observer & D50\n");
483 		spec = 1;
484 		illum = icxIT_D50;
485 		observ = icxOT_CIE_1931_2;
486 	}
487 
488 	/* Figure out what sort of device colorspace it is */
489 	{
490 		int ti;
491 
492 		if ((ti = icg->find_kword(icg, 0, "COLOR_REP")) < 0)
493 			error("Input file '%s' doesn't contain keyword COLOR_REPS",ti3name);
494 
495 		if (strcmp(icg->t[0].kdata[ti],"CMYK_XYZ") == 0) {
496 			devspace = icSigCmykData;
497 			devchan = 4;
498 			isLab = 0;
499 			isAdditive = 0;
500 		} else if (strcmp(icg->t[0].kdata[ti],"CMYK_LAB") == 0) {
501 			devspace = icSigCmykData;
502 			devchan = 4;
503 			isLab = 1;
504 			isAdditive = 0;
505 		} else if (strcmp(icg->t[0].kdata[ti],"CMY_XYZ") == 0) {
506 			devspace = icSigCmyData;
507 			devchan = 3;
508 			isLab = 0;
509 			isAdditive = 0;
510 		} else if (strcmp(icg->t[0].kdata[ti],"CMY_LAB") == 0) {
511 			devspace = icSigCmyData;
512 			devchan = 3;
513 			isLab = 1;
514 			isAdditive = 0;
515 		} else if (strcmp(icg->t[0].kdata[ti],"RGB_XYZ") == 0) {
516 			devspace = icSigRgbData;
517 			devchan = 3;
518 			isLab = 0;
519 			isAdditive = 1;
520 		} else if (strcmp(icg->t[0].kdata[ti],"RGB_LAB") == 0) {
521 			devspace = icSigRgbData;
522 			devchan = 3;
523 			isLab = 1;
524 			isAdditive = 1;
525 		} else if (strcmp(icg->t[0].kdata[ti],"iRGB_XYZ") == 0) {
526 			devspace = icSigRgbData;
527 			devchan = 3;
528 			isLab = 0;
529 			isAdditive = 1;
530 		} else if (strcmp(icg->t[0].kdata[ti],"iRGB_LAB") == 0) {
531 			devspace = icSigRgbData;
532 			devchan = 3;
533 			isLab = 1;
534 			isAdditive = 1;
535 		/* Scanner .ti3 files: */
536 		} else if (strcmp(icg->t[0].kdata[ti],"XYZ_RGB") == 0) {
537 			devspace = icSigRgbData;
538 			devchan = 3;
539 			isLab = 0;
540 			isAdditive = 1;
541 		} else if (strcmp(icg->t[0].kdata[ti],"LAB_RGB") == 0) {
542 			devspace = icSigRgbData;
543 			devchan = 3;
544 			isLab = 1;
545 			isAdditive = 1;
546 #ifdef IMP_MONO
547 		} else if (strcmp(icg->t[0].kdata[ti],"K_XYZ") == 0) {
548 			devspace = icSigGrayData;
549 			devchan = 1;
550 			isLab = 0;
551 			isAdditive = 0;
552 		} else if (strcmp(icg->t[0].kdata[ti],"K_LAB") == 0) {
553 			devspace = icSigGrayData;
554 			devchan = 1;
555 			isLab = 1;
556 			isAdditive = 0;
557 		} else if (strcmp(icg->t[0].kdata[ti],"W_XYZ") == 0) {
558 			devspace = icSigGrayData;
559 			devchan = 1;
560 			isLab = 0;
561 			isAdditive = 1;
562 		} else if (strcmp(icg->t[0].kdata[ti],"W_LAB") == 0) {
563 			devspace = icSigGrayData;
564 			devchan = 1;
565 			isLab = 1;
566 			isAdditive = 1;
567 #endif /* IMP_MONO */
568 
569 		} else
570 			error("Device input file '%s' has unhandled color representation '%s'",
571 			                                                     ti3name, icg->t[0].kdata[ti]);
572 	}
573 
574 	if ((npat = icg->t[0].nsets) <= 0)
575 		error("Input file '%s' has no sets of data",ti3name);
576 
577 	if (verb) {
578 		fprintf(verbo,"No of test patches = %d\n",npat);
579 	}
580 
581 	/* Allocate arrays to hold test patch input and output values */
582 	if ((tpat = (pval *)malloc(sizeof(pval) * npat)) == NULL)
583 		error("Malloc failed - tpat[]");
584 
585 	if ((stpat = (pval **)malloc(sizeof(pval *) * npat)) == NULL)
586 		error("Malloc failed - stpat[]");
587 
588 	for (i = 0; i < npat; i++)
589 		stpat[i] = &tpat[i];
590 
591 	/* Read in the CGATs fields */
592 	{
593 		int sidx;					/* Sample ID index */
594 		int sloc;					/* Sample location index (if any) */
595 		int ti, ci, mi, yi, ki;
596 		int Xi, Yi, Zi;
597 
598 		if ((sidx = icg->find_field(icg, 0, "SAMPLE_ID")) < 0)
599 			error("Input file '%s' doesn't contain field SAMPLE_ID",ti3name);
600 		if (icg->t[0].ftype[sidx] != nqcs_t)
601 			error("Input file '%s' field SAMPLE_ID is wrong type",ti3name);
602 
603 		if ((sloc = icg->find_field(icg, 0, "SAMPLE_LOC")) >= 0) {
604 			if (icg->t[0].ftype[sloc] != cs_t)
605 				error("Input file '%s' field SAMPLE_LOC is wrong type",ti3name);
606 		}
607 
608 		if (devspace == icSigGrayData) {
609 			if (isAdditive) {
610 				if ((ci = icg->find_field(icg, 0, "GRAY_W")) < 0)
611 					error("Input file doesn't contain field GRAY_W");
612 				if (icg->t[0].ftype[ci] != r_t)
613 					error("Field GRAY_W is wrong type - expect float");
614 			} else {
615 				if ((ci = icg->find_field(icg, 0, "GRAY_K")) < 0)
616 					error("Input file doesn't contain field GRAY_K");
617 				if (icg->t[0].ftype[ci] != r_t)
618 					error("Field GRAY_K is wrong type - expect float");
619 			}
620 			mi = yi = ki = ci;
621 
622 		} else if (devspace == icSigRgbData) {
623 			if ((ci = icg->find_field(icg, 0, "RGB_R")) < 0)
624 				error("Input file '%s' doesn't contain field RGB_R",ti3name);
625 			if (icg->t[0].ftype[ci] != r_t)
626 				error("Input file '%s' field RGB_R is wrong type - expect float",ti3name);
627 			if ((mi = icg->find_field(icg, 0, "RGB_G")) < 0)
628 				error("Input file '%s' doesn't contain field RGB_G",ti3name);
629 			if (icg->t[0].ftype[mi] != r_t)
630 				error("Input file '%s' field RGB_G is wrong type - expect float",ti3name);
631 			if ((yi = icg->find_field(icg, 0, "RGB_B")) < 0)
632 				error("Input file '%s' doesn't contain field RGB_B",ti3name);
633 			if (icg->t[0].ftype[yi] != r_t)
634 				error("Input file '%s' field RGB_B is wrong type - expect float",ti3name);
635 			ki = yi;
636 
637 		} else if (devspace == icSigCmyData) {
638 
639 			if ((ci = icg->find_field(icg, 0, "CMY_C")) < 0)
640 				error("Input file '%s' doesn't contain field CMY_C",ti3name);
641 			if (icg->t[0].ftype[ci] != r_t)
642 				error("Input file '%s' field CMY_C is wrong type - expect float",ti3name);
643 			if ((mi = icg->find_field(icg, 0, "CMY_M")) < 0)
644 				error("Input file '%s' doesn't contain field CMY_M",ti3name);
645 			if (icg->t[0].ftype[mi] != r_t)
646 				error("Input file '%s' field CMY_M is wrong type - expect float",ti3name);
647 			if ((yi = icg->find_field(icg, 0, "CMY_Y")) < 0)
648 				error("Input file '%s' doesn't contain field CMY_Y",ti3name);
649 			ki = yi;
650 		} else {	/* Assume CMYK */
651 
652 			if ((ci = icg->find_field(icg, 0, "CMYK_C")) < 0)
653 				error("Input file '%s' doesn't contain field CMYK_C",ti3name);
654 			if (icg->t[0].ftype[ci] != r_t)
655 				error("Input file '%s' field CMYK_C is wrong type - expect float",ti3name);
656 			if ((mi = icg->find_field(icg, 0, "CMYK_M")) < 0)
657 				error("Input file '%s' doesn't contain field CMYK_M",ti3name);
658 			if (icg->t[0].ftype[mi] != r_t)
659 				error("Input file '%s' field CMYK_M is wrong type - expect float",ti3name);
660 			if ((yi = icg->find_field(icg, 0, "CMYK_Y")) < 0)
661 				error("Input file '%s' doesn't contain field CMYK_Y",ti3name);
662 			if (icg->t[0].ftype[yi] != r_t)
663 				error("Input file '%s' field CMYK_Y is wrong type - expect float",ti3name);
664 			if ((ki = icg->find_field(icg, 0, "CMYK_K")) < 0)
665 				error("Input file '%s' doesn't contain field CMYK_K",ti3name);
666 			if (icg->t[0].ftype[ki] != r_t)
667 				error("Input file '%s' field CMYK_K is wrong type - expect float",ti3name);
668 		}
669 
670 		if (spec == 0) { 		/* Using instrument tristimulous value */
671 
672 			if (isLab) {		/* Expect Lab */
673 				if ((Xi = icg->find_field(icg, 0, "LAB_L")) < 0)
674 					error("Input file '%s' doesn't contain field LAB_L",ti3name);
675 				if (icg->t[0].ftype[Xi] != r_t)
676 					error("Input file '%s' field LAB_L is wrong type - expect float",ti3name);
677 				if ((Yi = icg->find_field(icg, 0, "LAB_A")) < 0)
678 					error("Input '%s' file doesn't contain field LAB_A",ti3name);
679 				if (icg->t[0].ftype[Yi] != r_t)
680 					error("Input file '%s' field LAB_A is wrong type - expect float",ti3name);
681 				if ((Zi = icg->find_field(icg, 0, "LAB_B")) < 0)
682 					error("Input file '%s' doesn't contain field LAB_B",ti3name);
683 				if (icg->t[0].ftype[Zi] != r_t)
684 					error("Input file '%s' field LAB_B is wrong type - expect float",ti3name);
685 
686 			} else { 		/* Expect XYZ */
687 				if ((Xi = icg->find_field(icg, 0, "XYZ_X")) < 0)
688 					error("Input file '%s' doesn't contain field XYZ_X",ti3name);
689 				if (icg->t[0].ftype[Xi] != r_t)
690 					error("Input file '%s' field XYZ_X is wrong type - expect float",ti3name);
691 				if ((Yi = icg->find_field(icg, 0, "XYZ_Y")) < 0)
692 					error("Input file '%s' doesn't contain field XYZ_Y",ti3name);
693 				if (icg->t[0].ftype[Yi] != r_t)
694 					error("Input file '%s' field XYZ_Y is wrong type - expect float",ti3name);
695 				if ((Zi = icg->find_field(icg, 0, "XYZ_Z")) < 0)
696 					error("Input file '%s' doesn't contain field XYZ_Z",ti3name);
697 				if (icg->t[0].ftype[Zi] != r_t)
698 					error("Input file '%s' field XYZ_Z is wrong type - expect float",ti3name);
699 			}
700 
701 			for (i = 0; i < npat; i++) {
702 				strcpy(tpat[i].sid, (char *)icg->t[0].fdata[i][sidx]);
703 				if (sloc >= 0)
704 					strcpy(tpat[i].slo, (char *)icg->t[0].fdata[i][sloc]);
705 				else
706 					strcpy(tpat[i].slo, "");
707 				tpat[i].p[0] = *((double *)icg->t[0].fdata[i][ci]) / 100.0;
708 				tpat[i].p[1] = *((double *)icg->t[0].fdata[i][mi]) / 100.0;
709 				tpat[i].p[2] = *((double *)icg->t[0].fdata[i][yi]) / 100.0;
710 				tpat[i].p[3] = *((double *)icg->t[0].fdata[i][ki]) / 100.0;
711 				if (tpat[i].p[0] > 1.0
712 				 || tpat[i].p[1] > 1.0
713 				 || tpat[i].p[2] > 1.0
714 				 || tpat[i].p[3] > 1.0) {
715 					error("Input file '%s' device value field value exceeds 100.0 !",ti3name);
716 				}
717 				tpat[i].v[0] = *((double *)icg->t[0].fdata[i][Xi]);
718 				tpat[i].v[1] = *((double *)icg->t[0].fdata[i][Yi]);
719 				tpat[i].v[2] = *((double *)icg->t[0].fdata[i][Zi]);
720 				if (!isLab && (!isdisp || isdnormed != 0)) {
721 					tpat[i].v[0] /= 100.0;		/* Normalise XYZ to range 0.0 - 1.0 */
722 					tpat[i].v[1] /= 100.0;
723 					tpat[i].v[2] /= 100.0;
724 				}
725 				if (!isLab) { /* Convert test patch result XYZ to PCS (D50 Lab) */
726 					icmXYZ2Lab(&icmD50, tpat[i].v, tpat[i].v);
727 				}
728 			}
729 
730 		} else { 		/* Using spectral data */
731 			int ii;
732 			xspect sp;
733 			char buf[100];
734 			int  spi[XSPECT_MAX_BANDS];	/* CGATS indexes for each wavelength */
735 			xsp2cie *sp2cie;	/* Spectral conversion object */
736 
737 			if ((ii = icg->find_kword(icg, 0, "SPECTRAL_BANDS")) < 0)
738 				error ("Input file '%s' doesn't contain keyword SPECTRAL_BANDS",ti3name);
739 			sp.spec_n = atoi(icg->t[0].kdata[ii]);
740 			if ((ii = icg->find_kword(icg, 0, "SPECTRAL_START_NM")) < 0)
741 				error ("Input file '%s' doesn't contain keyword SPECTRAL_START_NM",ti3name);
742 			sp.spec_wl_short = atof(icg->t[0].kdata[ii]);
743 			if ((ii = icg->find_kword(icg, 0, "SPECTRAL_END_NM")) < 0)
744 				error ("Input file '%s; doesn't contain keyword SPECTRAL_END_NM",ti3name);
745 			sp.spec_wl_long = atof(icg->t[0].kdata[ii]);
746 			if (!isdisp || isdnormed != 0)
747 				sp.norm = 100.0;
748 			else
749 				sp.norm = 1.0;
750 
751 			/* Find the fields for spectral values */
752 			for (j = 0; j < sp.spec_n; j++) {
753 				int nm;
754 
755 				/* Compute nearest integer wavelength */
756 				nm = (int)(sp.spec_wl_short + ((double)j/(sp.spec_n-1.0))
757 				            * (sp.spec_wl_long - sp.spec_wl_short) + 0.5);
758 
759 				sprintf(buf,"SPEC_%03d",nm);
760 
761 				if ((spi[j] = icg->find_field(icg, 0, buf)) < 0)
762 					error("Input file '%s' doesn't contain field %s",ti3name,buf);
763 
764 				if (icg->t[0].ftype[spi[j]] != r_t)
765 					error("Field %s is wrong type - expect float",buf);
766 			}
767 
768 			if (isemis) {
769 				illum = icxIT_none;		/* Emissive data has no illuminant */
770 				fwacomp = 0;
771 			}
772 
773 			/* Create a spectral conversion object */
774 			if ((sp2cie = new_xsp2cie(illum, illum == icxIT_none ? NULL : &cust_illum,
775 			                          observ, NULL, icSigLabData, icxClamp)) == NULL)
776 				error("Creation of spectral conversion object failed");
777 
778 			if (fwacomp) {
779 				double nw = 0.0;		/* Number of media white patches */
780 				xspect mwsp;			/* Medium spectrum */
781 				instType itype;			/* Spectral instrument type */
782 				xspect insp;			/* Instrument illuminant */
783 
784 				mwsp = sp;		/* Struct copy */
785 
786 	 			if ((ti = icg->find_kword(icg, 0, "TARGET_INSTRUMENT")) < 0)
787 					error ("Input file '%s' can't find target instrument needed for FWA compensation",ti3name);
788 
789 				if ((itype = inst_enum(icg->t[0].kdata[ti])) == instUnknown)
790 					error ("Input file '%s' unrecognised target instrument '%s'",ti3name, icg->t[0].kdata[ti]);
791 
792 				if (inst_illuminant(&insp, itype) != 0)
793 					error ("Instrument doesn't have an FWA illuminent");
794 
795 				/* Find the media white spectral reflectance */
796 				for (j = 0; j < mwsp.spec_n; j++)
797 					mwsp.spec[j] = 0.0;
798 
799 				/* Compute the mean of all the media white patches */
800 				for (i = 0; i < npat; i++) {
801 					int use = 0;
802 
803 					if (devspace == icSigGrayData) {
804 						if (isAdditive) {
805 							if (*((double *)icg->t[0].fdata[i][ci]) > (100.0 - 0.1))
806 								use = 1;
807 						} else {
808 							if (*((double *)icg->t[0].fdata[i][ci]) < 0.1)
809 								use = 1;
810 						}
811 					} else if (devspace == icSigRgbData) {
812 						if (*((double *)icg->t[0].fdata[i][ci]) > (100.0 - 0.1)
813 						 && *((double *)icg->t[0].fdata[i][mi]) > (100.0 - 0.1)
814 						 && *((double *)icg->t[0].fdata[i][yi]) > (100.0 - 0.1))
815 							use = 1;
816 					} else if (devspace == icSigCmyData) {
817 						if (*((double *)icg->t[0].fdata[i][ci]) < 0.1
818 						 && *((double *)icg->t[0].fdata[i][mi]) < 0.1
819 						 && *((double *)icg->t[0].fdata[i][yi]) < 0.1)
820 							use = 1;
821 					} else {	/* Assume CMYK */
822 
823 						if (*((double *)icg->t[0].fdata[i][ci]) < 0.1
824 						 && *((double *)icg->t[0].fdata[i][mi]) < 0.1
825 						 && *((double *)icg->t[0].fdata[i][yi]) < 0.1
826 						 && *((double *)icg->t[0].fdata[i][ki]) < 0.1) {
827 							use = 1;
828 						}
829 					}
830 
831 					if (use) {
832 						/* Read the spectral values for this patch */
833 						for (j = 0; j < mwsp.spec_n; j++) {
834 							mwsp.spec[j] += *((double *)icg->t[0].fdata[i][spi[j]]);
835 						}
836 						nw++;
837 					}
838 				}
839 				if (nw == 0.0) {
840 					warning("Input file '%s' can't find a media white patch to init FWA",ti3name);
841 
842 					/* Track the maximum reflectance for any band to determine white. */
843 					/* This might give bogus results if there is no white patch... */
844 					for (i = 0; i < npat; i++) {
845 						for (j = 0; j < mwsp.spec_n; j++) {
846 							double rv = *((double *)icg->t[0].fdata[i][spi[j]]);
847 							if (rv > mwsp.spec[j])
848 								mwsp.spec[j] = rv;
849 						}
850 					}
851 					nw++;
852 				}
853 
854 				for (j = 0; j < mwsp.spec_n; j++)
855 					mwsp.spec[j] /= nw;	/* Compute average */
856 
857 				/* If we are setting a specific simulated instrument illuminant */
858 				if (tillum != icxIT_none) {
859 					tillump = &cust_tillum;
860 					if (tillum != icxIT_custom) {
861 						if (standardIlluminant(tillump, tillum, 0.0)) {
862 							error("simulated inst. illum. not recognised");
863 						}
864 					}
865 				}
866 
867 				if (sp2cie->set_fwa(sp2cie, &insp, tillump, &mwsp))
868 					error ("Set FWA on sp2cie failed");
869 
870 				if (verb) {
871 					double FWAc;
872 					sp2cie->get_fwa_info(sp2cie, &FWAc);
873 					fprintf(verbo,"FWA content = %f\n",FWAc);
874 				}
875 
876 			}
877 
878 			for (i = 0; i < npat; i++) {
879 				strcpy(tpat[i].sid, (char *)icg->t[0].fdata[i][sidx]);
880 				if (sloc >= 0)
881 					strcpy(tpat[i].slo, (char *)icg->t[0].fdata[i][sloc]);
882 				else
883 					strcpy(tpat[i].slo, "");
884 				tpat[i].p[0] = *((double *)icg->t[0].fdata[i][ci]) / 100.0;
885 				tpat[i].p[1] = *((double *)icg->t[0].fdata[i][mi]) / 100.0;
886 				tpat[i].p[2] = *((double *)icg->t[0].fdata[i][yi]) / 100.0;
887 				tpat[i].p[3] = *((double *)icg->t[0].fdata[i][ki]) / 100.0;
888 				if (tpat[i].p[0] > 1.0
889 				 || tpat[i].p[1] > 1.0
890 				 || tpat[i].p[2] > 1.0
891 				 || tpat[i].p[3] > 1.0) {
892 					error("Input file '%s' device value field value exceeds 100.0 !",ti3name);
893 				}
894 
895 				/* Read the spectral values for this patch */
896 				for (j = 0; j < sp.spec_n; j++) {
897 					sp.spec[j] = *((double *)icg->t[0].fdata[i][spi[j]]);
898 				}
899 
900 				/* Convert it to CIE space */
901 				sp2cie->convert(sp2cie, tpat[i].v, &sp);
902 			}
903 
904 			sp2cie->del(sp2cie);		/* Done with this */
905 		}
906 		/* Normalize display values to Y = 1.0 if needed */
907 		/* (re-norm spec derived, since observer may be different) */
908 		if (isdisp && (isdnormed == 0 || spec != 0)) {
909 			double scale = -1e6;
910 			double bxyz[3];
911 
912 			/* Locate max Y */
913 			for (i = 0; i < npat; i++) {
914 				icmLab2XYZ(&icmD50, bxyz,  tpat[i].v);
915 				if (bxyz[1] > scale)
916 					scale = bxyz[1];
917 			}
918 
919 			scale = 1.0/scale;
920 
921 			/* Scale max Y to 1.0 */
922 			for (i = 0; i < npat; i++) {
923 				icmLab2XYZ(&icmD50, tpat[i].v, tpat[i].v);
924 				tpat[i].v[0] *= scale;
925 				tpat[i].v[1] *= scale;
926 				tpat[i].v[2] *= scale;
927 				icmXYZ2Lab(&icmD50, tpat[i].v, tpat[i].v);
928 			}
929 		}
930 
931 	}	/* End of reading in CGATs file */
932 
933 	/* - - - - - - - - - - */
934 	/* Compute the delta E of each point against the profile value */
935 
936 	/* Open up the file for reading */
937 	if ((rd_fp = new_icmFileStd_name(iccname,"r")) == NULL)
938 		error("Write: Can't open file '%s'",iccname);
939 
940 	if ((rd_icco = new_icc()) == NULL)
941 		error("Read: Creation of ICC object failed");
942 
943 	/* Read the header and tag list */
944 	if ((rv = rd_icco->read(rd_icco,rd_fp,0)) != 0)
945 		error("Read: %d, %s",rv,rd_icco->err);
946 
947 	/* Get the Fwd table, absolute with Lab override */
948 	if ((luo = rd_icco->get_luobj(rd_icco, icmFwd, intent,
949 	                              icSigLabData, icmLuOrdNorm)) == NULL) {
950 		error("%d, %s",rd_icco->errc, rd_icco->err);
951 	}
952 
953 	for (i = 0; i < npat; i++) {
954 
955 		/* Lookup the patch value in the profile */
956 		if (luo->lookup(luo, tpat[i].pv, tpat[i].p) > 1)
957 			error("%d, %s",rd_icco->errc,rd_icco->err);
958 
959 		if (cie2k)
960 			tpat[i].de = icmCIE2K(tpat[i].v, tpat[i].pv);
961 		else if (cie94)
962 			tpat[i].de = icmCIE94(tpat[i].v, tpat[i].pv);
963 		else
964 			tpat[i].de = icmLabDE(tpat[i].v, tpat[i].pv);
965 	}
966 
967 	/* - - - - - - - - - - */
968 	/* Sort by delta E */
969 	if (sortbyde) {
970 		/* Sort the dE's */
971 #define HEAP_COMPARE(A,B) (A.de > B.de)
972 		HEAPSORT(pval, tpat, npat);
973 #undef HEAP_COMPARE
974 	}
975 
976 	/* - - - - - - - - - - */
977 	/* Plot a dE histogram */
978 	if (dohisto) {
979 		double demax = -1e6, demin = 1e6;
980 		int maxbins = 50;		/* Maximum bins */
981 		int minbins = 20;		/* Target minimum bins (depends on distribution) */
982 		int mincount = 10;		/* Minimum number of points in a bin */
983 		double mbwidth;
984 		int nbins = 0;
985 		hbin *bins = NULL;
986 		double tval;
987 		double *x, *y;
988 
989 		/* Figure out the range of dE's */
990 		for (i = 0; i < npat; i++) {
991 			double de = tpat[i].de;
992 
993 			if (de > demax)
994 				demax = de;
995 			if (de < demin)
996 				demin = de;
997 		}
998 
999 		if (demax < 1e-6)
1000 			error("histogram: dE range is too small to plot");
1001 
1002 		/* Bin width that gives maxbins */
1003 		mbwidth = demax / maxbins;
1004 
1005 		/* Reduce mincount if needed to get minbins */
1006 		if (npat/minbins < mincount)
1007 			mincount = npat/minbins;
1008 
1009 		if ((bins = (hbin *)calloc(maxbins, sizeof(hbin))) == NULL)
1010 			error("malloc of histogram bins failed");
1011 
1012 		/* Sort the dE's */
1013 #define HEAP_COMPARE(A,B) (A->de < B->de)
1014 		HEAPSORT(pval *, stpat, npat);
1015 #undef HEAP_COMPARE
1016 
1017 		/* Create bins and add points */
1018 		bins[0].min = 0.0;
1019 		for (nbins = i = 0; i < npat && nbins < maxbins; i++) {
1020 			double de = stpat[i]->de;
1021 
1022 			/* Move on to next bin ? */
1023 			if (bins[nbins].count >= mincount
1024 			 && (de - bins[nbins].min) >= mbwidth) {
1025 				if (i > 0)
1026 					bins[nbins].max = 0.5 * (de + stpat[i-1]->de);
1027 				else
1028 					bins[nbins].max = de;
1029 				nbins++;
1030 				bins[nbins].min = bins[nbins-1].max;
1031 			}
1032 			bins[nbins].count++;
1033 			bins[nbins].max = de;
1034 		}
1035 		if (bins[nbins].count != 0)
1036 			nbins++;
1037 
1038 		/* Compute value */
1039 		tval = 0.0;
1040 		for (i = 0; i < nbins; i++) {
1041 			bins[i].val = bins[i].count/(bins[i].max - bins[i].min);
1042 			tval += bins[i].val;
1043 		}
1044 
1045 		tval /= 100.0;		/* Make it % */
1046 		for (i = 0; i < nbins; i++) {
1047 			bins[i].val /= tval;
1048 			if (verb) fprintf(verbo,"Bin %d, %f - %f, % 2.4f%%, count %d\n",
1049 			             i,bins[i].min,bins[i].max,bins[i].val,bins[i].count);
1050 		}
1051 
1052 		/* Plot it */
1053 		if ((x = (double *)calloc(nbins+1, sizeof(double))) == NULL)
1054 			error("malloc of histogram x array");
1055 		if ((y = (double *)calloc(nbins+1, sizeof(double))) == NULL)
1056 			error("malloc of histogram y array");
1057 
1058 		for (i = 0; i < nbins; i++) {
1059 			x[i] = 0.5 * (bins[i].min + bins[i].max);
1060 			y[i] = bins[i].val;
1061 		}
1062 		x[i] = demax;
1063 		y[i] = 0.0;
1064 		do_plot(x, y, NULL, NULL, nbins+1);
1065 
1066 		free(y);
1067 		free(x);
1068 		free(bins);
1069 	}
1070 
1071 	/* - - - - - - - - - - */
1072 	/* Create a pruned .ti3 file */
1073 	if (prune > 0.0) {
1074 		char *cp, outname[MAXNAMEL+31];
1075 		cgats *ocg;
1076 		cgats_set_elem *setel;		/* Array of set value elements */
1077 
1078 		strcpy(outname, ti3name);
1079 		if ((cp = strrchr(outname, '.')) == NULL)
1080 			cp = outname + strlen(outname);
1081 		sprintf(cp, "_p%.2f.ti3",prune);
1082 
1083 		if (verb) fprintf(verbo,"Created pruned file '%s'\n",outname);
1084 
1085 		/* Create the output files */
1086 		if ((ocg = new_cgats()) == NULL)
1087 			error("Failed to create cgats object");
1088 
1089 		/* Duplicate the type of the file */
1090 		if (icg->t[0].tt == cgats_X) {
1091 			ocg->add_other(ocg, icg->cgats_type);
1092 			ocg->add_table(ocg, tt_other, 0);
1093 		} else if (icg->t[0].tt == tt_other) {
1094 			ocg->add_other(ocg, icg->others[icg->t[0].oi]);
1095 			ocg->add_table(ocg, tt_other, 0);
1096 		} else {
1097 			ocg->add_table(ocg, icg->t[0].tt, 0);
1098 		}
1099 
1100 		/* Duplicate all the keywords */
1101 		for (i = 0; i < icg->t[0].nkwords; i++) {
1102 			ocg->add_kword(ocg, 0, icg->t[0].ksym[i], icg->t[0].kdata[i], NULL);
1103 		}
1104 
1105 		/* Duplicate all of the fields */
1106 		for (i = 0; i < icg->t[0].nfields; i++) {
1107 			ocg->add_field(ocg, 0, icg->t[0].fsym[i], icg->t[0].ftype[i]);
1108 		}
1109 
1110 		if ((setel = (cgats_set_elem *)malloc(
1111 		     sizeof(cgats_set_elem) * icg->t[0].nfields)) == NULL)
1112 			error("Malloc failed!");
1113 
1114 		/* Copy them approproately */
1115 		for (i = 0; i < icg->t[0].nsets; i++) {
1116 
1117 			if (tpat[i].de <= prune) {
1118 				icg->get_setarr(icg, 0, i, setel);
1119 				ocg->add_setarr(ocg, 0, setel);
1120 			}
1121 		}
1122 
1123 		if (verb) {
1124 			double acc;
1125 
1126 			acc = (double)ocg->t[0].nsets/(double)icg->t[0].nsets * 100.0;
1127 			fprintf(verbo,"%.2f%% accepted, %.3f%% rejected\n",acc, 100.0-acc);
1128 		}
1129 
1130 		/* Write out the file */
1131 		if (ocg->write_name(ocg, outname))
1132 			error("CGATS file '%s' write error : %s",outname,ocg->err);
1133 	}
1134 
1135 	/* - - - - - - - - - - */
1136 	/* Display various results */
1137 	{
1138 		double merr = 0.0;		/* Max */
1139 		double aerr = 0.0;		/* Avg */
1140 		double rerr = 0.0;		/* RMS */
1141 		double nsamps = 0.0;
1142 		int inn, outn;			/* Chanells for input and output spaces */
1143 
1144 		if (dovrml) {
1145 			wrl = new_vrml(out_name, doaxes, vrml_lab);
1146 			wrl->start_line_set(wrl, 0);
1147 		}
1148 
1149 		/* Get details of conversion (Arguments may be NULL if info not needed) */
1150 		luo->spaces(luo, NULL, &inn, NULL, &outn, NULL, NULL, NULL, NULL, NULL);
1151 
1152 		for (i = 0; i < npat; i++) {
1153 			double de, *out;
1154 
1155 			de = tpat[i].de;
1156 			out = tpat[i].pv;
1157 
1158 			if (verb > 1) {
1159 #ifdef HACK		// Print XYZ
1160 #pragma message("!!!!!!!!!!!!!!!!!! profile/profcheck.c HACK enabled !!!!!!!!!!!!!!!")
1161 				double outxyz[3], vxyz[3];
1162 				icmLab2XYZ(&icmD50, outxyz, out);
1163 				icmLab2XYZ(&icmD50, vxyz, tpat[i].v);
1164 
1165 				printf("[%f] %s%s%s: %s -> %f %f %f should be %f %f %f\n",
1166 				       de,
1167 				       tpat[i].sid,
1168 				       tpat[i].slo[0] != '\000' ? " @ " : "",
1169 				       tpat[i].slo,
1170 				       icmPdv(devchan, tpat[i].p),
1171 				       outxyz[0],outxyz[1],outxyz[2],
1172 				       vxyz[0],vxyz[1],vxyz[2]);
1173 #else
1174 				printf("[%f] %s%s%s: %s -> %f %f %f should be %f %f %f\n",
1175 				       de,
1176 				       tpat[i].sid,
1177 				       tpat[i].slo[0] != '\000' ? " @ " : "",
1178 				       tpat[i].slo,
1179 				       icmPdv(devchan, tpat[i].p),
1180 				       out[0],out[1],out[2],
1181 				       tpat[i].v[0],tpat[i].v[1],tpat[i].v[2]);
1182 #endif
1183 			}
1184 			if (dovrml) {
1185 				int ix[2];
1186 				double p1[3], p2[3];
1187 
1188 				/* Add the vertexes */
1189 				if (dominl && de < 0.5) {		/* Make a minimum length */
1190 					double cent[3], vec[3], vlen;
1191 
1192 					/* Compute center */
1193 					icmAdd3(cent, tpat[i].v, out);
1194 					icmScale3(cent, cent, 0.5);
1195 					if ((vlen = de) < 1e-6) {
1196 						vec[0] = 0.25; vec[1] = 0.0; vec[2] = 0.0;
1197 					} else {
1198 						icmSub3(vec, tpat[i].v, out);
1199 						icmScale3(vec, vec, 0.25/vlen);
1200 					}
1201 					icmSub3(p1, cent, vec);
1202 					icmAdd3(p2, cent, vec);
1203 				} else {
1204 					icmCpy3(p1, tpat[i].v);
1205 					icmCpy3(p2, out);
1206 				}
1207 				ix[0] = wrl->add_vertex(wrl, 0, p1);
1208 				ix[1] = wrl->add_vertex(wrl, 0, p2);
1209 
1210 				/* Add the line */
1211 				if (dodecol) {		/* Lines with color determined by length */
1212 					double rgb[3];
1213 					DE2RGB(rgb, icmNorm33(p1, p2));
1214 					wrl->add_col_line(wrl, 0, ix, rgb);
1215 
1216 				} else {	/* Natural color */
1217 					wrl->add_line(wrl, 0, ix);
1218 				}
1219 			}
1220 
1221 			/* Stats */
1222 			aerr += de;
1223 			rerr += de * de;
1224 
1225 			nsamps++;
1226 			if (de > merr)
1227 				merr = de;
1228 
1229 		}
1230 		if (dovrml) {
1231 			wrl->make_lines_vc(wrl, 0, 0.0);
1232 			wrl->del(wrl);
1233 		}
1234 		printf("Profile check complete, errors%s: max. = %f, avg. = %f, RMS = %f\n",
1235             cie2k ? "(CIEDE2000)" : cie94 ? " (CIE94)" : "", merr, aerr/nsamps, sqrt(rerr/nsamps));
1236 
1237 		/* ------------------------------- */
1238 		/* If we want sort by target value */
1239 		if (ddevv) {
1240 			double cieval[3];
1241 
1242 			/* Lookup the CIE value of the target */
1243 			if (luo->lookup(luo, cieval, devval) > 1)
1244 				error("%d, %s",rd_icco->errc,rd_icco->err);
1245 
1246 			/* Compute deltas to target value. */
1247 			for (i = 0; i < npat; i++) {
1248 				if (cie2k)
1249 					tpat[i].dv = icmCIE2K(tpat[i].v, cieval);
1250 				else if (cie94)
1251 					tpat[i].dv = icmCIE94(tpat[i].v, cieval);
1252 				else
1253 					tpat[i].dv = icmLabDE(tpat[i].v, cieval);
1254 
1255 				tpat[i].dp = 0.0;
1256 				for (j = 0; j < inn; j++) {
1257 					double tt;
1258 					tt = tpat[i].p[j] - devval[j];
1259 					tpat[i].dp += tt * tt;
1260 				}
1261 				tpat[i].dp = sqrt(tpat[i].dp);
1262 			}
1263 
1264 			if (sortbypcs) {
1265 				/* Sort by pcs delta */
1266 #define HEAP_COMPARE(A,B) (A->dv < B->dv)
1267 				HEAPSORT(pval *, stpat, npat);
1268 #undef HEAP_COMPARE
1269 			} else {
1270 				/* Sort by device delta */
1271 #define HEAP_COMPARE(A,B) (A->dp < B->dp)
1272 				HEAPSORT(pval *, stpat, npat);
1273 #undef HEAP_COMPARE
1274 			}
1275 
1276 			printf("Target point:\n");
1277 			if (devspace == icSigCmykData) {
1278 				printf("%f %f %f %f -> %f %f %f\n",
1279 				       devval[0],devval[1],devval[2],devval[3],
1280 				       cieval[0],cieval[1],cieval[2]);
1281 			} else {	/* Assume RGB/CMY */
1282 				printf("%f %f %f -> %f %f %f\n",
1283 				       devval[0],devval[1],devval[2],
1284 				       cieval[0],cieval[1],cieval[2]);
1285 			}
1286 			printf("\n");
1287 
1288 			for (i = 0; i < npat; i++) {
1289 				if (devspace == icSigCmykData) {
1290 					printf("%s: %f %f %f %f [%f] -> %f %f %f [%f]\n",
1291 					       stpat[i]->sid,
1292 					       stpat[i]->p[0],stpat[i]->p[1],stpat[i]->p[2],stpat[i]->p[3],
1293 					       stpat[i]->dp,
1294 					       stpat[i]->v[0],stpat[i]->v[1],stpat[i]->v[2],
1295 					       stpat[i]->dv);
1296 				} else {	/* Assume RGB/CMY */
1297 					printf("%s: %f %f %f [%f] -> %f %f %f [%f]\n",
1298 					       stpat[i]->sid,
1299 					       stpat[i]->p[0],stpat[i]->p[1],stpat[i]->p[2],
1300 					       stpat[i]->dp,
1301 					       stpat[i]->v[0],stpat[i]->v[1],stpat[i]->v[2],
1302 					       stpat[i]->dv);
1303 				}
1304 			}
1305 		}
1306 
1307 		/* Done with lookup object */
1308 		luo->del(luo);
1309 
1310 		/* Close the file */
1311 		rd_icco->del(rd_icco);
1312 		rd_fp->del(rd_fp);
1313 	}
1314 
1315 	icg->del(icg);		/* Clean up */
1316 
1317 	return 0;
1318 }
1319 
1320 
1321 /* ------------------------------------------------ */
1322 
1323 /* Convert a delta E value into a signal color: */
1324 static void
DE2RGB(double * out,double in)1325 DE2RGB(double *out, double in) {
1326 	struct {
1327 		double de;
1328 		double r, g, b;
1329 	} range[6] = {
1330 		{ 10.0, 1, 1, 0 },		/* yellow */
1331 		{ 4.0,  1, 0, 0 },		/* red */
1332 		{ 2.0, 1, 0, 1 },		/* magenta */
1333 		{ 1.0, 0, 0, 1 },		/* blue */
1334 		{ 0.5, 0, 1, 1 },		/* cyan */
1335 		{ 0.0, 0, 1, 0 }		/* green */
1336 	};
1337 	int i;
1338 	double bl;
1339 
1340 	/* Locate the range we're in */
1341 	if (in > range[0].de) {
1342 		out[0] = range[0].r;
1343 		out[1] = range[0].g;
1344 		out[2] = range[0].b;
1345 	} else {
1346 		for (i = 0; i < 5; i++) {
1347 			if (in <= range[i].de && in >= range[i+1].de)
1348 				break;
1349 		}
1350 		bl = (in - range[i+1].de)/(range[i].de - range[i+1].de);
1351 		out[0] = bl * range[i].r + (1.0 - bl) * range[i+1].r;
1352 		out[1] = bl * range[i].g + (1.0 - bl) * range[i+1].g;
1353 		out[2] = bl * range[i].b + (1.0 - bl) * range[i+1].b;
1354 	}
1355 }
1356 
1357 
1358 
1359