1470 lines
44 KiB
C
1470 lines
44 KiB
C
/*
|
|
* Argyll Color Management System
|
|
* Verify two sets of PCS values.
|
|
*
|
|
* Author: Graeme W. Gill
|
|
* Date: 7/6/2005
|
|
*
|
|
* Copyright 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 two CGATS files (probably but not necesserily .ti3 files) of PCS
|
|
* values (either XYZ, L*a*b* or spectral), matches the values, and computes
|
|
* overall errors. This is useful for verifying proofing systems.
|
|
*/
|
|
|
|
/*
|
|
* TTBD:
|
|
*
|
|
* We're no warning about setting illuminant or FWA for emissive spectral data,
|
|
* they are just ignored.
|
|
*/
|
|
|
|
#undef DEBUG
|
|
|
|
#define verbo stdout
|
|
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
#include "copyright.h"
|
|
#include "aconfig.h"
|
|
#include "numlib.h"
|
|
#include "vrml.h"
|
|
#include "cgats.h"
|
|
#include "conv.h"
|
|
#include "xicc.h"
|
|
#include "insttypes.h"
|
|
#include "disptechs.h"
|
|
#include "ccmx.h"
|
|
#include "sort.h"
|
|
#include "plot.h"
|
|
#include "ui.h"
|
|
|
|
#ifdef DEBUG
|
|
#undef DBG
|
|
#define DBG(xxx) printf xxx ;
|
|
#else
|
|
#undef DBG
|
|
#define DBG(xxx)
|
|
#endif
|
|
|
|
void
|
|
usage(void) {
|
|
fprintf(stderr,"Verify CIE values, Version %s\n",ARGYLL_VERSION_STR);
|
|
fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
|
|
fprintf(stderr,"usage: colverify [-options] target.ti3 measured.ti3\n");
|
|
fprintf(stderr," -v [n] Verbose mode, n >= 2 print each value\n");
|
|
fprintf(stderr," -l Match patches by sample location rather than id\n");
|
|
fprintf(stderr," -n Normalise each files reading to its white Y\n");
|
|
fprintf(stderr," -N Normalise each files reading to its white XYZ\n");
|
|
fprintf(stderr," -m Normalise each files reading to its white X+Y+Z\n");
|
|
fprintf(stderr," -M Normalise both files reading to mean white XYZ\n");
|
|
fprintf(stderr," -D Use D50 100.0 as L*a*b* white reference\n");
|
|
fprintf(stderr," -c Show CIE94 delta E values\n");
|
|
fprintf(stderr," -k Show CIEDE2000 delta E values\n");
|
|
fprintf(stderr," -h [hist.txt] Plot a histogram of delta E's [Optionally save points to .txt]\n");
|
|
fprintf(stderr," -s Sort patch values by error\n");
|
|
fprintf(stderr," -w create PCS %s vector visualisation (measured%s)\n",vrml_format(),vrml_ext());
|
|
fprintf(stderr," -W create PCS %s marker visualisation (measured%s)\n",vrml_format(),vrml_ext());
|
|
fprintf(stderr," -d create Device RGB %s marker visualisation (measured%s)\n",vrml_format(),vrml_ext());
|
|
// fprintf(stderr," -d y create Device YCbCr %s marker visualisation (measured%s)\n",vrml_format(),vrml_ext());
|
|
fprintf(stderr," -x Use %s axes\n",vrml_format());
|
|
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, 2015_2, 2015_10, S&B 1955_2, shaw, J&V 1978_2 or file.cmf\n");
|
|
fprintf(stderr," -L profile.%s Skip any first file, out of profile gamut patches\n",ICC_FILE_EXT_ND);
|
|
fprintf(stderr," -X file.ccmx Apply Colorimeter Correction Matrix to second file\n");
|
|
// fprintf(stderr," -Z A|X Just print Average|Max +tab\n");
|
|
fprintf(stderr," target.ti3 Target (reference) PCS or spectral values.\n");
|
|
fprintf(stderr," measured.ti3 Measured (actual) PCS or spectral values\n");
|
|
exit(1);
|
|
}
|
|
|
|
/* Patch value type */
|
|
typedef struct {
|
|
char sid[50]; /* sample id */
|
|
char loc[100]; /* sample location (empty if none) */
|
|
double rgb[3]; /* RGB value if RGB device space present, or YCbCr if dovrml==4 */
|
|
double ycc[3]; /* YCbCr if RGB and dovrml==4 */
|
|
int og; /* Out of gamut flag */
|
|
double xyz[3]; /* XYZ value */
|
|
double v[3]; /* Lab value */
|
|
double de; /* Delta E */
|
|
double ixde[3]; /* XYZ Component DE */
|
|
double ide[3]; /* Lab Component DE */
|
|
} pval;
|
|
|
|
/* Histogram bin type */
|
|
typedef struct {
|
|
int count; /* Raw count */
|
|
double val; /* Normalized value */
|
|
double min, max; /* Bin range */
|
|
} hbin;
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
int fa,nfa,mfa; /* current argument we're looking at */
|
|
int verb = 0; /* Verbose level */
|
|
int useloc = 0; /* Match patches by sample location */
|
|
int norm = 0; /* 1 = norm to White Y, 2 = norm to White XYZ */
|
|
/* 3 = norm to White X+Y+Z, 4 = norm to average XYZ */
|
|
int usestdd50 = 0; /* Use standard D50 instead of avg white as reference */
|
|
int cie94 = 0;
|
|
int cie2k = 0;
|
|
int dovrml = 0; /* 1 = PCS vector, 2 = PCS marker, 3 = RGB, 4 - YCbCr */
|
|
int doaxes = 0;
|
|
int dohisto = 0; /* Plot histogram of delta E's */
|
|
static char histoname[MAXNAMEL+1] = "\000"; /* Optional file to save histogram points to */
|
|
int dosort = 0;
|
|
int dozrep = 0; /* 1 = print average, 2 = print max */
|
|
static char ccmxname[MAXNAMEL+1] = "\000"; /* Colorimeter Correction Matrix name */
|
|
ccmx *cmx = NULL; /* Colorimeter Correction Matrix */
|
|
static char gprofname[MAXNAMEL+1] = "\000"; /* Gamut limit profile name */
|
|
icmFile *fp = NULL;
|
|
icmErr err = { 0, { '\000'} };
|
|
icc *icco = NULL;
|
|
xicc *xicco = NULL;
|
|
icxLuBase *luo = NULL;
|
|
|
|
struct {
|
|
char name[MAXNAMEL+1]; /* Patch filename */
|
|
int isemis; /* nz if emsissive spectral reference data */
|
|
int isdisp; /* nz if display */
|
|
int isdnormed; /* Has display data been normalised to 100 ? */
|
|
int isrgb; /* Is RGB device space ? */
|
|
int npat; /* Number of patches */
|
|
int nig; /* Number of patches in gamut */
|
|
double w[3]; /* XYZ of "white" */
|
|
double nw[3]; /* Normalised XYZ of "white" */
|
|
pval *pat; /* patch values */
|
|
} cg[2]; /* Target and current patch file information */
|
|
|
|
int *match; /* Array mapping first list indexes to corresponding second */
|
|
int *sort; /* Array of first list indexes in sorted order */
|
|
int fwacomp = 0; /* FWA compensation */
|
|
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_none; /* Spectral defaults to D50 */
|
|
xspect cust_illum; /* Custom illumination spectrum */
|
|
icxObserverType obType = icxOT_none; /* Defaults to 1931 2 degree */
|
|
xspect custObserver[3]; /* If obType = icxOT_custom */
|
|
|
|
icmXYZNumber labw = icmD50; /* The Lab white reference */
|
|
|
|
static char out_name[MAXNAMEL+4+1]; /* VRML/X3D name */
|
|
vrml *wrl = NULL;
|
|
|
|
|
|
int i, j, n;
|
|
|
|
if (argc <= 1)
|
|
usage();
|
|
|
|
/* Process the arguments */
|
|
mfa = 2; /* Minimum final 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+mfa) < argc) {
|
|
if (argv[fa+1][0] != '-') {
|
|
nfa = fa + 1;
|
|
na = argv[nfa]; /* next is seperate non-flag argument */
|
|
}
|
|
}
|
|
}
|
|
|
|
if (argv[fa][1] == '?')
|
|
usage();
|
|
|
|
/* Verbose */
|
|
else if (argv[fa][1] == 'v') {
|
|
verb = 1;
|
|
|
|
if (na != NULL && na[0] >= '0' && na[0] <= '9') {
|
|
fa = nfa;
|
|
verb = atoi(na);
|
|
}
|
|
}
|
|
|
|
/* Use location to match patches */
|
|
else if (argv[fa][1] == 'l') {
|
|
useloc = 1;
|
|
}
|
|
|
|
/* normalize */
|
|
else if (argv[fa][1] == 'n'
|
|
|| argv[fa][1] == 'N') {
|
|
norm = 1;
|
|
if (argv[fa][1] == 'N')
|
|
norm = 2;
|
|
}
|
|
|
|
else if (argv[fa][1] == 'm') {
|
|
norm = 3;
|
|
}
|
|
|
|
else if (argv[fa][1] == 'M') {
|
|
norm = 4;
|
|
}
|
|
|
|
else if (argv[fa][1] == 'D')
|
|
usestdd50 = 1;
|
|
|
|
/* VRML/X3D */
|
|
else if (argv[fa][1] == 'w')
|
|
dovrml = 1;
|
|
|
|
else if (argv[fa][1] == 'W')
|
|
dovrml = 2;
|
|
|
|
else if (argv[fa][1] == 'd') {
|
|
dovrml = 3;
|
|
if (na != NULL) { /* Argument is present - RGB or YCbCr. */
|
|
fa = nfa;
|
|
if (strcmp(na, "y") == 0)
|
|
dovrml = 4;
|
|
else
|
|
usage();
|
|
}
|
|
}
|
|
|
|
/* Axes */
|
|
else if (argv[fa][1] == 'x')
|
|
doaxes = 1;
|
|
|
|
/* CIE94 delta E */
|
|
else if (argv[fa][1] == 'c') {
|
|
cie94 = 1;
|
|
cie2k = 0;
|
|
}
|
|
|
|
else if (argv[fa][1] == 'k') {
|
|
cie94 = 0;
|
|
cie2k = 1;
|
|
}
|
|
|
|
/* Plot histogram */
|
|
else if (argv[fa][1] == 'h') {
|
|
dohisto = 1;
|
|
|
|
if (na != NULL) { /* Argument is present - file to save points to */
|
|
fa = nfa;
|
|
strncpy(histoname,na,MAXNAMEL); histoname[MAXNAMEL] = '\000';
|
|
}
|
|
}
|
|
|
|
/* Sort */
|
|
else if (argv[fa][1] == 's')
|
|
dosort = 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 */
|
|
inst_meas_type mt;
|
|
|
|
spec = 1;
|
|
tillum = icxIT_custom;
|
|
if (read_xspect(&cust_tillum, &mt, NULL, na) != 0)
|
|
usage();
|
|
|
|
if (mt != inst_mrt_none
|
|
&& mt != inst_mrt_emission
|
|
&& mt != inst_mrt_ambient
|
|
&& mt != inst_mrt_emission_flash
|
|
&& mt != inst_mrt_ambient_flash) {
|
|
error("Target illuminant '%s' is wrong measurement type",na);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Spectral to CIE Illuminant type */
|
|
else if (argv[fa][1] == 'i') {
|
|
if (na == NULL) usage();
|
|
fa = nfa;
|
|
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 */
|
|
inst_meas_type mt;
|
|
|
|
spec = 1;
|
|
illum = icxIT_custom;
|
|
if (read_xspect(&cust_illum, &mt, NULL, na) != 0)
|
|
usage();
|
|
|
|
if (mt != inst_mrt_none
|
|
&& mt != inst_mrt_emission
|
|
&& mt != inst_mrt_ambient
|
|
&& mt != inst_mrt_emission_flash
|
|
&& mt != inst_mrt_ambient_flash) {
|
|
error("CIE illuminant '%s' is wrong measurement type",na);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Spectral Observer type */
|
|
else if (argv[fa][1] == 'o') {
|
|
if (na == NULL) usage();
|
|
fa = nfa;
|
|
if (strcmp(na, "1931_2") == 0) { /* Classic 2 degree */
|
|
spec = 1;
|
|
obType = icxOT_CIE_1931_2;
|
|
} else if (strcmp(na, "1964_10") == 0) { /* Classic 10 degree */
|
|
spec = 1;
|
|
obType = icxOT_CIE_1964_10;
|
|
} else if (strcmp(na, "2015_2") == 0) { /* Latest 2 degree */
|
|
spec = 1;
|
|
obType = icxOT_CIE_2015_2;
|
|
} else if (strcmp(na, "2015_10") == 0) { /* Latest 10 degree */
|
|
spec = 1;
|
|
obType = icxOT_CIE_2015_10;
|
|
} else if (strcmp(na, "1955_2") == 0) { /* Stiles and Burch 1955 2 degree */
|
|
spec = 1;
|
|
obType = icxOT_Stiles_Burch_2;
|
|
} else if (strcmp(na, "1978_2") == 0) { /* Judd and Voss 1978 2 degree */
|
|
spec = 1;
|
|
obType = icxOT_Judd_Voss_2;
|
|
} else if (strcmp(na, "shaw") == 0) { /* Shaw and Fairchilds 1997 2 degree */
|
|
spec = 1;
|
|
obType = icxOT_Shaw_Fairchild_2;
|
|
} else { /* Assume it's a filename */
|
|
obType = icxOT_custom;
|
|
if (read_cmf(custObserver, na) != 0)
|
|
usage();
|
|
}
|
|
}
|
|
|
|
/* Gamut limit profile for first file */
|
|
else if (argv[fa][1] == 'L') {
|
|
if (na == NULL) usage();
|
|
fa = nfa;
|
|
strncpy(gprofname,na,MAXNAMEL-1); gprofname[MAXNAMEL-1] = '\000';
|
|
}
|
|
|
|
/* Colorimeter Correction Matrix for second file */
|
|
else if (argv[fa][1] == 'X') {
|
|
if (na == NULL) usage();
|
|
fa = nfa;
|
|
strncpy(ccmxname,na,MAXNAMEL-1); ccmxname[MAXNAMEL-1] = '\000';
|
|
}
|
|
|
|
else if (argv[fa][1] == 'Z') {
|
|
if (na == NULL) usage();
|
|
fa = nfa;
|
|
if (strcmp(na, "A") == 0) {
|
|
dozrep = 1;
|
|
} else if (strcmp(na, "X") == 0) {
|
|
dozrep = 2;
|
|
} else
|
|
usage();
|
|
|
|
} else
|
|
usage();
|
|
} else
|
|
break;
|
|
}
|
|
|
|
/* Get the file name arguments */
|
|
if (fa >= argc || argv[fa][0] == '-') usage();
|
|
strncpy(cg[0].name,argv[fa++],MAXNAMEL); cg[0].name[MAXNAMEL] = '\000';
|
|
|
|
if (fa >= argc || argv[fa][0] == '-') usage();
|
|
strncpy(cg[1].name,argv[fa],MAXNAMEL); cg[1].name[MAXNAMEL] = '\000';
|
|
|
|
|
|
/* Create VRML/X3D base name */
|
|
{
|
|
char *xl;
|
|
strcpy(out_name, cg[1].name);
|
|
if ((xl = strrchr(out_name, '.')) == NULL) /* Figure where extention is */
|
|
xl = out_name + strlen(out_name);
|
|
xl[0] = '\000'; /* Remove extension */
|
|
}
|
|
|
|
if (fwacomp && spec == 0)
|
|
error ("FWA compensation only works when viewer and/or illuminant selected");
|
|
|
|
/* Gamut limit profile */
|
|
if (gprofname[0] != '\000') {
|
|
int rv;
|
|
|
|
if ((fp = new_icmFileStd_name(&err,gprofname,"r")) == NULL)
|
|
error ("Can't open file '%s' (0x%x, '%s')",gprofname,err.c,err.m);
|
|
|
|
if ((icco = new_icc(&err)) == NULL)
|
|
error ("Creation of ICC object failed (0x%x, '%s')",err.c,err.m);
|
|
|
|
if ((rv = icco->read(icco,fp,0)) != 0)
|
|
error("Reading profile '%s' failed failed with error %d:'%s'\n",
|
|
gprofname, icco->e.c, icco->e.m);
|
|
|
|
if (icco->header->deviceClass != icSigInputClass
|
|
&& icco->header->deviceClass != icSigDisplayClass
|
|
&& icco->header->deviceClass != icSigOutputClass)
|
|
error("Profile '%s' must be a device profile to filter by gamut",gprofname);
|
|
|
|
/* Wrap with an expanded icc */
|
|
if ((xicco = new_xicc(icco)) == NULL)
|
|
error ("Creation of xicc failed");
|
|
|
|
/* Get a expanded color conversion object */
|
|
if ((luo = xicco->get_luobj(xicco, ICX_CLIP_NEAREST | ICX_FAST_SETUP,
|
|
icmFwd, icRelativeColorimetric, icSigXYZData, icmLuOrdNorm, NULL, NULL)) == NULL)
|
|
error ("%d, %s",xicco->e.c, xicco->e.m);
|
|
}
|
|
|
|
/* Colorimeter Correction Matrix */
|
|
if (ccmxname[0] != '\000') {
|
|
if ((cmx = new_ccmx()) == NULL)
|
|
error("new_ccmx failed\n");
|
|
if (cmx->read_ccmx(cmx,ccmxname))
|
|
error("Reading Colorimeter Correction Matrix file '%s' failed with error %d:'%s'\n",
|
|
ccmxname, cmx->e.c, cmx->e.m);
|
|
}
|
|
|
|
|
|
/* Open up each file in turn, target then measured, */
|
|
/* and read in the CIE values. */
|
|
for (n = 0; n < 2; n++) {
|
|
cgats *cgf = NULL; /* cgats file data */
|
|
int isLab = 0; /* 0 if file CIE is XYZ, 1 if is Lab */
|
|
int sidx; /* Sample ID index */
|
|
int sldx = -1; /* Sample location index, < 0 if invalid */
|
|
int xix, yix, zix;
|
|
int rgbix[3]; /* RGB field indexes (if rgb ) */
|
|
int dti; /* Device type index */
|
|
|
|
/* Open CIE target values */
|
|
cgf = new_cgats(); /* Create a CGATS structure */
|
|
cgf->add_other(cgf, ""); /* Allow any signature file */
|
|
|
|
DBG(("Opening file '%s'\n",cg[n].name))
|
|
|
|
if (cgf->read_name(cgf, cg[n].name))
|
|
error("CGATS file '%s' read error : %s",cg[n].name,cgf->e.m);
|
|
|
|
if (cgf->ntables < 1)
|
|
error ("Input file '%s' doesn't contain at least one table",cg[n].name);
|
|
|
|
if ((dti = cgf->find_kword(cgf, 0, "DEVICE_CLASS")) < 0)
|
|
warning("Input file '%s' doesn't contain keyword DEVICE_CLASS",cg[n].name);
|
|
|
|
/* Figure out what sort of device it is */
|
|
{
|
|
int ti;
|
|
|
|
cg[n].isemis = 0;
|
|
cg[n].isdisp = 0;
|
|
cg[n].isdnormed = 0;
|
|
cg[n].w[0] = cg[n].w[1] = cg[n].w[2] = 0.0;
|
|
|
|
if (dti >= 0) {
|
|
if (strcmp(cgf->t[0].kdata[dti],"DISPLAY") == 0) {
|
|
cg[n].isemis = 1;
|
|
cg[n].isdisp = 1;
|
|
cg[n].isdnormed = 1; /* Assume display type is normalised to 100 */
|
|
|
|
} else if (strcmp(cgf->t[0].kdata[dti],"EMISINPUT") == 0) {
|
|
cg[n].isemis = 1;
|
|
}
|
|
|
|
if (cg[n].isdisp) {
|
|
|
|
if ((ti = cgf->find_kword(cgf, 0, "LUMINANCE_XYZ_CDM2")) >= 0) {
|
|
if (sscanf(cgf->t[0].kdata[ti], " %lf %lf %lf ",&cg[n].w[0], &cg[n].w[1], &cg[n].w[2]) != 3)
|
|
cg[n].w[0] = cg[n].w[1] = cg[n].w[2] = 0.0;
|
|
}
|
|
|
|
/* See if there is an explicit tag indicating data has been normalised to Y = 100 */
|
|
if ((ti = cgf->find_kword(cgf, 0, "NORMALIZED_TO_Y_100")) >= 0) {
|
|
if (strcmp(cgf->t[0].kdata[ti],"NO") == 0) {
|
|
cg[n].isdnormed = 0;
|
|
} else {
|
|
cg[n].isdnormed = 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Check if the file is suitable */
|
|
if (!spec
|
|
&& cgf->find_field(cgf, 0, "LAB_L") < 0
|
|
&& cgf->find_field(cgf, 0, "XYZ_X") < 0) {
|
|
|
|
if (cgf->find_kword(cgf, 0, "SPECTRAL_BANDS") < 0)
|
|
error ("Neither CIE nor spectral data found in file '%s'",cg[n].name);
|
|
|
|
/* Switch to using spectral information */
|
|
if (verb)
|
|
printf("No CIE data found, switching to spectral with standard observer & D50 for file '%s'\n",cg[n].name);
|
|
spec = 1;
|
|
}
|
|
if (spec && cgf->find_kword(cgf, 0, "SPECTRAL_BANDS") < 0)
|
|
error ("No spectral data data found in file '%s' when spectral expected",cg[n].name);
|
|
|
|
if (!spec && cgf->find_field(cgf, 0, "LAB_L") >= 0)
|
|
isLab = 1;
|
|
|
|
cg[n].nig = cg[n].npat = cgf->t[0].nsets; /* Number of patches */
|
|
|
|
/* See if it has RGB device space (for -d option) */
|
|
if ((rgbix[0] = cgf->find_field(cgf, 0, "RGB_R")) >= 0
|
|
&& cgf->t[0].ftype[rgbix[0]] == r_t
|
|
|
|
&& (rgbix[1] = cgf->find_field(cgf, 0, "RGB_G")) >= 0
|
|
&& cgf->t[0].ftype[rgbix[1]] == r_t
|
|
|
|
&& (rgbix[2] = cgf->find_field(cgf, 0, "RGB_B")) >= 0
|
|
&& cgf->t[0].ftype[rgbix[2]] == r_t) {
|
|
cg[n].isrgb = 1;
|
|
} else {
|
|
cg[n].isrgb = 0;
|
|
}
|
|
|
|
/* Read all the target patches */
|
|
if (cg[n].npat <= 0)
|
|
error("No sets of data in file '%s'",cg[n].name);
|
|
|
|
if (verb && n == 0) {
|
|
fprintf(verbo,"No of test patches = %d\n",cg[n].npat);
|
|
}
|
|
|
|
/* Allocate arrays to hold test patch input and output values */
|
|
if ((cg[n].pat = (pval *)malloc(sizeof(pval) * cg[n].npat)) == NULL)
|
|
error("Malloc failed - pat[]");
|
|
|
|
/* Read in the CGATs fields */
|
|
if ((sidx = cgf->find_field(cgf, 0, "SAMPLE_ID")) < 0
|
|
&& (sidx = cgf->find_field(cgf, 0, "SampleName")) < 0
|
|
&& (sidx = cgf->find_field(cgf, 0, "Sample_Name")) < 0
|
|
&& (sidx = cgf->find_field(cgf, 0, "SAMPLE_NAME")) < 0)
|
|
error("Input file '%s' doesn't contain field SAMPLE_ID, SampleName, Sample_Name, SAMPLE_NAME",cg[n].name);
|
|
if (cgf->t[0].ftype[sidx] != nqcs_t
|
|
&& cgf->t[0].ftype[sidx] != cs_t)
|
|
error("Sample ID/Name field isn't a quoted or non quoted character string");
|
|
|
|
if ((sldx = cgf->find_field(cgf, 0, "SAMPLE_LOC")) < 0
|
|
|| cgf->t[0].ftype[sldx] != cs_t)
|
|
sldx = -1;
|
|
|
|
if (spec == 0) { /* Using instrument tristimulous value */
|
|
|
|
if (isLab) { /* Expect Lab */
|
|
if ((xix = cgf->find_field(cgf, 0, "LAB_L")) < 0)
|
|
error("Input file '%s' doesn't contain field LAB_L",cg[n].name);
|
|
if (cgf->t[0].ftype[xix] != r_t)
|
|
error("Field LAB_L is wrong type - expect float");
|
|
if ((yix = cgf->find_field(cgf, 0, "LAB_A")) < 0)
|
|
error("Input file '%s' doesn't contain field LAB_A",cg[n].name);
|
|
if (cgf->t[0].ftype[yix] != r_t)
|
|
error("Field LAB_A is wrong type - expect float");
|
|
if ((zix = cgf->find_field(cgf, 0, "LAB_B")) < 0)
|
|
error("Input file '%s' doesn't contain field LAB_B",cg[n].name);
|
|
if (cgf->t[0].ftype[zix] != r_t)
|
|
error("Field LAB_B is wrong type - expect float");
|
|
|
|
} else { /* Expect XYZ */
|
|
if ((xix = cgf->find_field(cgf, 0, "XYZ_X")) < 0)
|
|
error("Input file '%s' doesn't contain field XYZ_X",cg[n].name);
|
|
if (cgf->t[0].ftype[xix] != r_t)
|
|
error("Field XYZ_X is wrong type - expect float");
|
|
if ((yix = cgf->find_field(cgf, 0, "XYZ_Y")) < 0)
|
|
error("Input file '%s' doesn't contain field XYZ_Y",cg[n].name);
|
|
if (cgf->t[0].ftype[yix] != r_t)
|
|
error("Field XYZ_Y is wrong type - expect float");
|
|
if ((zix = cgf->find_field(cgf, 0, "XYZ_Z")) < 0)
|
|
error("Input file '%s' doesn't contain field XYZ_Z",cg[n].name);
|
|
if (cgf->t[0].ftype[zix] != r_t)
|
|
error("Field XYZ_Z is wrong type - expect float");
|
|
}
|
|
|
|
for (i = 0; i < cg[n].npat; i++) {
|
|
strcpy(cg[n].pat[i].sid, (char *)cgf->t[0].fdata[i][sidx]);
|
|
if (sldx >= 0)
|
|
strcpy(cg[n].pat[i].loc, (char *)cgf->t[0].fdata[i][sldx]);
|
|
else
|
|
cg[n].pat[i].loc[0] = '\000';
|
|
cg[n].pat[i].og = 0;
|
|
cg[n].pat[i].xyz[0] = *((double *)cgf->t[0].fdata[i][xix]);
|
|
cg[n].pat[i].xyz[1] = *((double *)cgf->t[0].fdata[i][yix]);
|
|
cg[n].pat[i].xyz[2] = *((double *)cgf->t[0].fdata[i][zix]);
|
|
|
|
if (isLab) { /* Convert Lab to XYZ 0..100% */
|
|
icmLab2XYZ(&icmD50_100, cg[n].pat[i].xyz, cg[n].pat[i].xyz);
|
|
}
|
|
//printf("~1 file %d patch %d = XYZ %f %f %f\n", n,i,cg[n].pat[i].xyz[0],cg[n].pat[i].xyz[1],cg[n].pat[i].xyz[2]);
|
|
|
|
/* restore normalised display values to absolute */
|
|
if (cg[n].isdnormed) {
|
|
if (cg[n].w[1] > 0.0) { // Found absoluute display white tag
|
|
cg[n].pat[i].xyz[0] *= cg[n].w[1]/100.0;
|
|
cg[n].pat[i].xyz[1] *= cg[n].w[1]/100.0;
|
|
cg[n].pat[i].xyz[2] *= cg[n].w[1]/100.0;
|
|
}
|
|
|
|
} else if (!cg[n].isdisp) {
|
|
/* If reflective or transmissive that are 0..100%, */
|
|
/* scale back to 0.. 1 */
|
|
cg[n].pat[i].xyz[0] /= 100.0; /* scale back to XYZ 1.0 */
|
|
cg[n].pat[i].xyz[1] /= 100.0;
|
|
cg[n].pat[i].xyz[2] /= 100.0;
|
|
}
|
|
|
|
/* Apply ccmx */
|
|
if (n == 1 && cmx != NULL) {
|
|
cmx->xform(cmx, cg[n].pat[i].xyz, cg[n].pat[i].xyz);
|
|
}
|
|
|
|
if ((dovrml == 3 || dovrml == 4) && cg[n].isrgb) {
|
|
cg[n].pat[i].rgb[0] = 0.01 * *((double *)cgf->t[0].fdata[i][rgbix[0]]);
|
|
cg[n].pat[i].rgb[1] = 0.01 * *((double *)cgf->t[0].fdata[i][rgbix[1]]);
|
|
cg[n].pat[i].rgb[2] = 0.01 * *((double *)cgf->t[0].fdata[i][rgbix[2]]);
|
|
|
|
if (dovrml == 4) {
|
|
icmRec709_RGBd_2_YPbPr(cg[n].pat[i].ycc, cg[n].pat[i].rgb);
|
|
}
|
|
}
|
|
}
|
|
|
|
} 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 */
|
|
|
|
/* Copies of global values: */
|
|
int l_fwacomp = fwacomp;
|
|
int l_spec = spec;
|
|
icxIllumeType l_tillum = tillum;
|
|
xspect *l_tillump = tillump;
|
|
icxIllumeType l_illum = illum;
|
|
icxObserverType l_observ = obType;
|
|
|
|
if ((ii = cgf->find_kword(cgf, 0, "SPECTRAL_BANDS")) < 0)
|
|
error ("Input file doesn't contain keyword SPECTRAL_BANDS");
|
|
sp.spec_n = atoi(cgf->t[0].kdata[ii]);
|
|
if ((ii = cgf->find_kword(cgf, 0, "SPECTRAL_START_NM")) < 0)
|
|
error ("Input file doesn't contain keyword SPECTRAL_START_NM");
|
|
sp.spec_wl_short = atof(cgf->t[0].kdata[ii]);
|
|
if ((ii = cgf->find_kword(cgf, 0, "SPECTRAL_END_NM")) < 0)
|
|
error ("Input file doesn't contain keyword SPECTRAL_END_NM");
|
|
sp.spec_wl_long = atof(cgf->t[0].kdata[ii]);
|
|
if (!cg[n].isdisp || cg[n].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] = cgf->find_field(cgf, 0, buf)) < 0)
|
|
error("Input file doesn't contain field %s",buf);
|
|
|
|
if (cgf->t[0].ftype[spi[j]] != r_t)
|
|
error("Field %s is wrong type - expect float",buf);
|
|
}
|
|
|
|
/* Create a spectral conversion object */
|
|
if (cg[n].isemis) {
|
|
if (l_illum != icxIT_none)
|
|
warning("-i illuminant ignored for emissive reference type");
|
|
if (l_fwacomp)
|
|
warning("-f FWA ignored for emissive reference type");
|
|
l_illum = icxIT_none; /* Make emissive conversion */
|
|
l_tillum = icxIT_none;
|
|
l_fwacomp = 0;
|
|
} else {
|
|
/* Set default */
|
|
if (l_illum == icxIT_none)
|
|
l_illum = icxIT_D50;
|
|
}
|
|
|
|
/* Set default */
|
|
if (l_observ == icxOT_none)
|
|
l_observ = icxOT_CIE_1931_2;
|
|
|
|
if ((sp2cie = new_xsp2cie(l_illum, 0.0, l_illum == icxIT_none ? NULL : &cust_illum,
|
|
l_observ, custObserver, icSigXYZData, icxClamp)) == NULL)
|
|
error("Creation of spectral conversion object failed");
|
|
|
|
if (l_fwacomp) {
|
|
int ti;
|
|
xspect mwsp; /* Medium spectrum */
|
|
instType itype; /* Spectral instrument type */
|
|
xspect insp; /* Instrument illuminant */
|
|
|
|
mwsp = sp; /* Struct copy */
|
|
|
|
if ((ti = cgf->find_kword(cgf, 0, "TARGET_INSTRUMENT")) < 0)
|
|
error ("Can't find target instrument in '%s' needed for FWA compensation",cg[n].name);
|
|
|
|
if ((itype = inst_enum(cgf->t[0].kdata[ti])) == instUnknown)
|
|
error ("Unrecognised target instrument '%s'", cgf->t[0].kdata[ti]);
|
|
|
|
if (inst_illuminant(&insp, itype) != 0)
|
|
error ("Instrument doesn't have an FWA illuminent");
|
|
|
|
/* Determine a media white spectral reflectance */
|
|
for (j = 0; j < mwsp.spec_n; j++)
|
|
mwsp.spec[j] = 0.0;
|
|
|
|
/* Since we don't want to assume that there are any associated device */
|
|
/* values present in each file, we can't use this as means of */
|
|
/* determining the media color. Use an alternative approach here, */
|
|
/* which may give slightly different results to profile. */
|
|
|
|
/* Track the maximum reflectance for any band to determine white. */
|
|
/* This might silently fail, if there isn't white in the sample set. */
|
|
for (i = 0; i < cg[0].npat; i++) {
|
|
for (j = 0; j < mwsp.spec_n; j++) {
|
|
double rv = *((double *)cgf->t[0].fdata[i][spi[j]]);
|
|
if (rv > mwsp.spec[j])
|
|
mwsp.spec[j] = rv;
|
|
}
|
|
}
|
|
|
|
/* If we are setting a specific simulated instrument illuminant */
|
|
if (l_tillum != icxIT_none) {
|
|
l_tillump = &cust_tillum;
|
|
if (l_tillum != icxIT_custom) {
|
|
if (standardIlluminant(l_tillump, l_tillum, 0.0)) {
|
|
error("simulated inst. illum. not recognised");
|
|
}
|
|
}
|
|
}
|
|
|
|
if (sp2cie->set_fwa(sp2cie, &insp, l_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 < cg[0].npat; i++) {
|
|
|
|
strcpy(cg[n].pat[i].sid, (char *)cgf->t[0].fdata[i][sidx]);
|
|
if (sldx >= 0)
|
|
strcpy(cg[n].pat[i].loc, (char *)cgf->t[0].fdata[i][sldx]);
|
|
else
|
|
cg[n].pat[i].loc[0] = '\000';
|
|
cg[n].pat[i].og = 0;
|
|
|
|
/* Read the spectral values for this patch */
|
|
for (j = 0; j < sp.spec_n; j++) {
|
|
sp.spec[j] = *((double *)cgf->t[0].fdata[i][spi[j]]);
|
|
}
|
|
|
|
/* Convert it to XYZ space */
|
|
sp2cie->convert(sp2cie, cg[n].pat[i].xyz, &sp);
|
|
|
|
/* restore normalised display values to absolute */
|
|
if (cg[n].isdnormed) {
|
|
if (cg[n].w[1] > 0.0) { // Found absoluute display white tag
|
|
cg[n].pat[i].xyz[0] *= cg[n].w[1];
|
|
cg[n].pat[i].xyz[1] *= cg[n].w[1];
|
|
cg[n].pat[i].xyz[2] *= cg[n].w[1];
|
|
}
|
|
|
|
}
|
|
|
|
/* Apply ccmx */
|
|
if (n == 1 && cmx != NULL) {
|
|
cmx->xform(cmx, cg[n].pat[i].xyz, cg[n].pat[i].xyz);
|
|
}
|
|
}
|
|
|
|
sp2cie->del(sp2cie); /* Done with this */
|
|
} /* End of reading in CGATs file spectral data */
|
|
|
|
|
|
/* Locate the patch with maximum Y, a possible white patch */
|
|
/* in case we need it latter. */
|
|
{
|
|
int ii;
|
|
|
|
if (cg[n].w[1] == 0.0) { /* No display white patch tag */
|
|
|
|
/* Locate patch with biggest Y, assume it is white... */
|
|
for (i = 0; i < cg[n].npat; i++) {
|
|
if (cg[n].pat[i].xyz[1] > cg[n].w[1]) {
|
|
icmCpy3(cg[n].w, cg[n].pat[i].xyz);
|
|
ii = i;
|
|
}
|
|
}
|
|
if (verb) printf("File %d Chose patch %d as white, XYZ %f %f %f\n",
|
|
n, ii+1,cg[n].w[0],cg[n].w[1],cg[n].w[2]);
|
|
} else {
|
|
if (verb) printf("File %d White is from display luminance ref. XYZ %f %f %f\n",
|
|
n, cg[n].w[0],cg[n].w[1],cg[n].w[2]);
|
|
}
|
|
icmCpy3(cg[n].nw, cg[n].w);
|
|
}
|
|
cgf->del(cgf); /* Clean up */
|
|
} /* Next file */
|
|
|
|
if (norm == 4) { /* Normalise to average of white XYZ of the two files */
|
|
icmBlend3(cg[0].w, cg[0].w, cg[1].w, 0.5);
|
|
icmCpy3(cg[1].w, cg[0].w);
|
|
// if (verb) printf("Average White XYZ %f %f %f\n",cg[0].w[0],cg[0].w[1],cg[0].w[2]);
|
|
}
|
|
|
|
/* For both files */
|
|
for (n = 0; n < 2; n++) {
|
|
|
|
/* Normalise this file to white = 1.0 or D50 */
|
|
if (norm) {
|
|
int ii;
|
|
|
|
double chmat[3][3]; /* Chromatic adapation matrix */
|
|
|
|
DBG(("Normalizng '%s' to white\n",cg[n].name))
|
|
|
|
if (norm == 2 || norm == 4) { /* Norm to white XYZ */
|
|
icmXYZNumber s_wp;
|
|
icmAry2XYZ(s_wp, cg[n].w);
|
|
icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, s_wp, chmat);
|
|
}
|
|
|
|
for (i = 0; i < cg[n].npat; i++) {
|
|
if (norm == 1) {
|
|
cg[n].pat[i].xyz[0] *= 100.0 / cg[n].w[1];
|
|
cg[n].pat[i].xyz[1] *= 100.0 / cg[n].w[1];
|
|
cg[n].pat[i].xyz[2] *= 100.0 / cg[n].w[1];
|
|
} else if (norm == 2 || norm == 4) {
|
|
icmMulBy3x3(cg[n].pat[i].xyz, chmat, cg[n].pat[i].xyz);
|
|
} else {
|
|
cg[n].pat[i].xyz[0] *= 100.0 / (cg[n].w[0] + cg[n].w[1] + cg[n].w[2]);
|
|
cg[n].pat[i].xyz[1] *= 100.0 / (cg[n].w[0] + cg[n].w[1] + cg[n].w[2]);
|
|
cg[n].pat[i].xyz[2] *= 100.0 / (cg[n].w[0] + cg[n].w[1] + cg[n].w[2]);
|
|
}
|
|
//printf("~1 file %d patch %d = norm XYZ %f %f %f\n", n,i,cg[n].pat[i].xyz[0],cg[n].pat[i].xyz[1],cg[n].pat[i].xyz[2]);
|
|
}
|
|
/* Compute normalised white too */
|
|
if (norm == 1) {
|
|
cg[n].nw[0] *= 100.0 / cg[n].w[1];
|
|
cg[n].nw[1] *= 100.0 / cg[n].w[1];
|
|
cg[n].nw[2] *= 100.0 / cg[n].w[1];
|
|
} else if (norm == 2 || norm == 4) {
|
|
icmMulBy3x3(cg[n].nw, chmat, cg[n].w);
|
|
} else {
|
|
cg[n].nw[0] *= 100.0 / (cg[n].w[0] + cg[n].w[1] + cg[n].w[2]);
|
|
cg[n].nw[1] *= 100.0 / (cg[n].w[0] + cg[n].w[1] + cg[n].w[2]);
|
|
cg[n].nw[2] *= 100.0 / (cg[n].w[0] + cg[n].w[1] + cg[n].w[2]);
|
|
}
|
|
//printf("~1 file %d norm white XYZ %f %f %f\n", n,cg[n].nw[0], cg[n].nw[1], cg[n].nw[2]);
|
|
}
|
|
} /* Next file */
|
|
|
|
if (cmx != NULL)
|
|
cmx->del(cmx);
|
|
cmx = NULL;
|
|
|
|
/* Check that the number of test patches matches */
|
|
if (cg[0].npat != cg[1].npat)
|
|
error("Number of patches between '%s' and '%s' doesn't match",cg[0].name,cg[1].name);
|
|
|
|
/* Create a list to map the second list of patches to the first */
|
|
if ((match = (int *)malloc(sizeof(int) * cg[0].npat)) == NULL)
|
|
error("Malloc failed - match[]");
|
|
|
|
/* Use location to match */
|
|
if (useloc) {
|
|
for (i = 0; i < cg[0].npat; i++) {
|
|
if (cg[0].pat[0].loc[0] == '\000'
|
|
|| cg[1].pat[0].loc[0] == '\000')
|
|
error("Need both files to have SAMPLE_LOC field to match on location");
|
|
|
|
for (j = 0; j < cg[1].npat; j++) {
|
|
if (strcmp(cg[0].pat[i].loc, cg[1].pat[j].loc) == 0)
|
|
break; /* Found it */
|
|
}
|
|
if (j < cg[1].npat) {
|
|
match[i] = j;
|
|
} else {
|
|
error("Failed to find matching patch to '%s'",cg[0].pat[i].loc);
|
|
}
|
|
}
|
|
|
|
/* Use id */
|
|
} else {
|
|
for (i = 0; i < cg[0].npat; i++) {
|
|
for (j = 0; j < cg[1].npat; j++) {
|
|
if (strcmp(cg[0].pat[i].sid, cg[1].pat[j].sid) == 0)
|
|
break; /* Found it */
|
|
}
|
|
if (j < cg[1].npat) {
|
|
match[i] = j;
|
|
} else {
|
|
error("Failed to find matching patch to '%s'",cg[0].pat[i].sid);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Figure out which patches to skip because they are out of gamut */
|
|
if (luo != NULL) {
|
|
double chmat[3][3]; /* Chromatic adapation matrix */
|
|
double out[MAX_CHAN], in[3], check[3];
|
|
icmXYZNumber s_wp;
|
|
int rv;
|
|
|
|
DBG(("Figuring out of gamut patches\n"))
|
|
|
|
/* Convert sample PCS to relative */
|
|
//printf(" cg[0].w %f %f %f\n", cg[0].w[0], cg[0].w[1], cg[0].w[2]);
|
|
icmAry2XYZ(s_wp, cg[0].w);
|
|
s_wp.X /= s_wp.Y; // Normalise the white to 1.0
|
|
s_wp.Y /= s_wp.Y; // so that matrix doesn't change magnitude
|
|
s_wp.Z /= s_wp.Y;
|
|
//printf(" s_wp %f %f %f\n", s_wp.X, s_wp.Y, s_wp.Z);
|
|
icmChromAdaptMatrix(ICM_CAM_BRADFORD, icmD50, s_wp, chmat);
|
|
//printf("~1 matrix = \n");
|
|
//printf(" %f %f %f\n", chmat[0][0], chmat[0][1], chmat[0][2]);
|
|
//printf(" %f %f %f\n", chmat[1][0], chmat[1][1], chmat[1][2]);
|
|
//printf(" %f %f %f\n", chmat[2][0], chmat[2][1], chmat[2][2]);
|
|
|
|
for (i = 0; i < cg[0].npat; i++) {
|
|
icmMulBy3x3(in, chmat, cg[0].pat[i].xyz);
|
|
|
|
//printf("~1 %d: xyz %f %f %f, rel %f %f %f\n", i+1, cg[0].pat[i].xyz[0], cg[0].pat[i].xyz[1], cg[0].pat[i].xyz[2], in[0], in[1], in[2]);
|
|
|
|
if ((rv = luo->inv_lookup(luo, out, in)) > 0 || 1) {
|
|
double de;
|
|
|
|
luo->lookup(luo, check, out);
|
|
de = icmXYZLabDE(&icmD50,check, in);
|
|
//printf("~1 %d: rv %d, de %f, check XYZ %f %f %f\n",i+1,rv, de, check[0],check[1],check[2]);
|
|
|
|
if (de >= 0.01) {
|
|
cg[0].pat[i].og = 1;
|
|
//printf("~1 Patch %d is out of gamut by DE %f RGB %f %f %f\n",i+1,de,out[0],out[1],out[2]);
|
|
if (verb >= 3)
|
|
printf("Patch %d is out of gamut by DE %f\n",i+1,de);
|
|
cg[0].nig--;
|
|
}
|
|
}
|
|
}
|
|
if (verb)
|
|
fprintf(verbo,"No of test patches in gamut = %d/%d\n",cg[0].nig,cg[0].npat);
|
|
}
|
|
|
|
if (cg[0].nig <= 0) {
|
|
if (verb)
|
|
fprintf(verbo,"No test patches in gamut - givig up\n");
|
|
return 0;
|
|
}
|
|
|
|
/* Adjust the Lab reference white to be the mean of the white of the two files */
|
|
if (norm != 0 && !usestdd50) {
|
|
labw.X = 0.5 * (cg[0].nw[0] + cg[1].nw[0]);
|
|
labw.Y = 0.5 * (cg[0].nw[1] + cg[1].nw[1]);
|
|
labw.Z = 0.5 * (cg[0].nw[2] + cg[1].nw[2]);
|
|
|
|
if (verb)
|
|
printf("L*a*b* white reference = XYZ %f %f %f\n",labw.X,labw.Y,labw.Z);
|
|
}
|
|
|
|
/* labw defaults to D50 */
|
|
|
|
/* Convert XYZ to Lab */
|
|
for (n = 0; n < 2; n++) {
|
|
for (i = 0; i < cg[n].npat; i++) {
|
|
icmXYZ2Lab(&labw, cg[n].pat[i].v, cg[n].pat[i].xyz);
|
|
}
|
|
}
|
|
|
|
/* Compute the delta E's */
|
|
DBG(("Computing the delta E's\n"))
|
|
for (i = 0; i < cg[0].npat; i++) {
|
|
|
|
cg[0].pat[i].ixde[0] = fabs(cg[0].pat[i].xyz[0] - cg[1].pat[match[i]].xyz[0]);
|
|
cg[0].pat[i].ixde[1] = fabs(cg[0].pat[i].xyz[1] - cg[1].pat[match[i]].xyz[1]);
|
|
cg[0].pat[i].ixde[2] = fabs(cg[0].pat[i].xyz[2] - cg[1].pat[match[i]].xyz[2]);
|
|
|
|
if (cie2k)
|
|
cg[0].pat[i].de = icmCIE2K(cg[0].pat[i].v, cg[1].pat[match[i]].v);
|
|
else if (cie94)
|
|
cg[0].pat[i].de = icmCIE94(cg[0].pat[i].v, cg[1].pat[match[i]].v);
|
|
else
|
|
cg[0].pat[i].de = icmLabDE(cg[0].pat[i].v, cg[1].pat[match[i]].v);
|
|
|
|
cg[0].pat[i].ide[0] = fabs(cg[0].pat[i].v[0] - cg[1].pat[match[i]].v[0]);
|
|
cg[0].pat[i].ide[1] = fabs(cg[0].pat[i].v[1] - cg[1].pat[match[i]].v[1]);
|
|
cg[0].pat[i].ide[2] = fabs(cg[0].pat[i].v[2] - cg[1].pat[match[i]].v[2]);
|
|
}
|
|
|
|
/* Create sorted list, from worst to best. */
|
|
if ((sort = (int *)malloc(sizeof(int) * cg[0].npat)) == NULL)
|
|
error("Malloc failed - sort[]");
|
|
for (i = 0; i < cg[0].npat; i++)
|
|
sort[i] = i;
|
|
|
|
#define HEAP_COMPARE(A,B) (cg[0].pat[A].de > cg[0].pat[B].de)
|
|
HEAPSORT(int, sort, cg[0].npat);
|
|
#undef HEAP_COMPARE
|
|
|
|
/* - - - - - - - - - - */
|
|
/* Plot a dE histogram */
|
|
if (dohisto) {
|
|
double demax = -1e6, demin = 1e6;
|
|
int maxbins = 50; /* Maximum bins */
|
|
// int minbins = 20; /* Target minimum bins (depends on distribution) */
|
|
// int mincount = 10; /* Minimum number of points in a bin */
|
|
int minbins = 10; /* Target minimum bins (depends on distribution) */
|
|
int mincount = 5; /* Minimum number of points in a bin */
|
|
double mbwidth;
|
|
int nbins = 0;
|
|
hbin *bins = NULL;
|
|
pval **stpat; /* Pointers to sorted cg[0].pat[] */
|
|
double tval;
|
|
double *x, *y;
|
|
|
|
DBG(("Plotting histogram\n"))
|
|
|
|
/* Figure out the range of dE's */
|
|
for (i = 0; i < cg[0].npat; i++) {
|
|
double de = cg[0].pat[i].de;
|
|
|
|
if (de > demax)
|
|
demax = de;
|
|
if (de < demin)
|
|
demin = de;
|
|
}
|
|
|
|
if (demax < 1e-6)
|
|
error("histogram: dE range is too small to plot");
|
|
|
|
/* Bin width that gives maxbins */
|
|
mbwidth = demax / maxbins;
|
|
|
|
#ifdef NEVER
|
|
/* Reduce mincount if needed to get minbins */
|
|
if (cg[0].npat/minbins < mincount)
|
|
mincount = cg[0].npat/minbins;
|
|
#endif
|
|
|
|
if ((bins = (hbin *)calloc(maxbins, sizeof(hbin))) == NULL)
|
|
error("malloc of histogram bins failed");
|
|
|
|
if ((stpat = (pval **)malloc(sizeof(pval *) * cg[0].npat)) == NULL)
|
|
error("Malloc failed - stpat[]");
|
|
|
|
for (i = 0; i < cg[0].npat; i++)
|
|
stpat[i] = &cg[0].pat[i];
|
|
|
|
/* Sort the dE's */
|
|
#define HEAP_COMPARE(A,B) (A->de < B->de)
|
|
HEAPSORT(pval *, stpat, cg[0].npat);
|
|
#undef HEAP_COMPARE
|
|
|
|
/* Create bins and add points */
|
|
bins[0].min = 0.0;
|
|
for (nbins = i = 0; i < cg[0].npat && nbins < maxbins; i++) {
|
|
double de = stpat[i]->de;
|
|
|
|
/* Move on to next bin ? */
|
|
if (bins[nbins].count >= mincount
|
|
&& (de - bins[nbins].min) >= mbwidth) {
|
|
if (i > 0)
|
|
bins[nbins].max = 0.5 * (de + stpat[i-1]->de);
|
|
else
|
|
bins[nbins].max = de;
|
|
nbins++;
|
|
bins[nbins].min = bins[nbins-1].max;
|
|
}
|
|
bins[nbins].count++;
|
|
bins[nbins].max = de;
|
|
}
|
|
if (bins[nbins].count != 0)
|
|
nbins++;
|
|
|
|
/* Compute value */
|
|
tval = 0.0;
|
|
for (i = 0; i < nbins; i++) {
|
|
bins[i].val = bins[i].count/(bins[i].max - bins[i].min);
|
|
tval += bins[i].val;
|
|
}
|
|
|
|
tval /= 100.0; /* Make it % */
|
|
for (i = 0; i < nbins; i++) {
|
|
bins[i].val /= tval;
|
|
if (verb) fprintf(verbo,"Bin %d, %f - %f, % 2.4f%%, count %d\n",
|
|
i,bins[i].min,bins[i].max,bins[i].val,bins[i].count);
|
|
}
|
|
|
|
/* Plot it */
|
|
if ((x = (double *)calloc(nbins+1, sizeof(double))) == NULL)
|
|
error("malloc of histogram x array");
|
|
if ((y = (double *)calloc(nbins+1, sizeof(double))) == NULL)
|
|
error("malloc of histogram y array");
|
|
|
|
for (i = 0; i < nbins; i++) {
|
|
x[i] = 0.5 * (bins[i].min + bins[i].max);
|
|
y[i] = bins[i].val;
|
|
}
|
|
x[i] = demax;
|
|
y[i] = 0.0;
|
|
do_plot(x, y, NULL, NULL, nbins+1);
|
|
|
|
/* Save points to a file ? */
|
|
if (histoname[0] != '\000') {
|
|
FILE *fp;
|
|
|
|
if ((fp = fopen(histoname, "w")) == NULL)
|
|
error("Opening '%s' for writing failed",histoname);
|
|
|
|
fprintf(fp, "%s\t%s\n\n",cg[0].name, cg[1].name);
|
|
|
|
for (i = 0; i < nbins; i++) {
|
|
fprintf(fp, "%f\t%f\n",0.5 * (bins[i].min + bins[i].max), bins[i].val);
|
|
}
|
|
|
|
if (fclose(fp) != 0)
|
|
error("Closing '%s' failed",histoname);
|
|
}
|
|
|
|
free(y);
|
|
free(x);
|
|
free(bins);
|
|
free(stpat);
|
|
}
|
|
|
|
/* - - - - - - - - - - */
|
|
/* Figure out the report */
|
|
{
|
|
double merr = 0.0, aerr = 0.0;
|
|
int n90;
|
|
double merr90 = 0.0, aerr90 = 0.0;
|
|
int n10;
|
|
double merr10 = 0.0, aerr10 = 0.0;
|
|
double rad;
|
|
double aierr[3] = { 0.0, 0.0, 0.0 };
|
|
double aixerr[3] = { 0.0, 0.0, 0.0 };
|
|
double red[3] = { 1.0, 0.2, 0.2 };
|
|
double green[3] = { 0.2, 1.0, 0.2 };
|
|
double min[3], max[3];
|
|
double col[3];
|
|
|
|
if (dovrml) {
|
|
double vol;
|
|
int k;
|
|
|
|
wrl = new_vrml(out_name, doaxes, (dovrml == 3 || dovrml == 4) ? vrml_rgb : vrml_lab);
|
|
wrl->start_line_set(wrl, 0);
|
|
|
|
for (j = 0; j < 3; j++) {
|
|
min[j] = 1e6;
|
|
max[j] = -1e6;
|
|
}
|
|
|
|
/* Get bounding box */
|
|
for (i = 0; i < cg[0].npat; i++) {
|
|
for (k = 0; k < 2; k++) {
|
|
for (j = 0; j < 3; j++) {
|
|
if (dovrml == 3 || dovrml == 4) { /* RGB or YCC device plot */
|
|
if (cg[k].pat[i].rgb[j] > max[j])
|
|
max[j] = cg[k].pat[i].rgb[j];
|
|
if (cg[k].pat[i].rgb[j] < min[j])
|
|
min[j] = cg[k].pat[i].rgb[j];
|
|
} else {
|
|
if (cg[k].pat[i].v[j] > max[j])
|
|
max[j] = cg[k].pat[i].v[j];
|
|
if (cg[k].pat[i].v[j] < min[j])
|
|
min[j] = cg[k].pat[i].v[j];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
for (vol = 1.0, j = 0; j < 3; j++) {
|
|
//printf("~1 size[%d] = %f\n",j, max[j] - min[j]);
|
|
vol *= (max[j] - min[j]);
|
|
}
|
|
vol = sqrt(vol);
|
|
//printf("~1 vol = %f\n",vol);
|
|
rad = 0.02 * vol/pow(cg[0].npat, 1.0/3.0);
|
|
//printf("~1 rad = %f\n",rad);
|
|
|
|
if (dovrml == 3) // Hack
|
|
rad = 0.02;
|
|
else if (dovrml == 4) // Hack
|
|
rad = 0.015;
|
|
}
|
|
|
|
if (dovrml && (dovrml == 3 || dovrml == 4)) { /* RGB/YCC device plot */
|
|
if (!cg[0].isrgb || !cg[1].isrgb)
|
|
error("Both files must have RGB devices space for -d option");
|
|
}
|
|
|
|
/* Do overall results */
|
|
for (i = 0; i < cg[0].npat; i++) {
|
|
double de;
|
|
|
|
if (cg[0].pat[i].og) /* Skip out of gamut patches */
|
|
continue;
|
|
|
|
if (dosort)
|
|
j = sort[i];
|
|
else
|
|
j = i;
|
|
|
|
de = cg[0].pat[j].de;
|
|
aerr += de;
|
|
|
|
aierr[0] += cg[0].pat[j].ide[0];
|
|
aierr[1] += cg[0].pat[j].ide[1];
|
|
aierr[2] += cg[0].pat[j].ide[2];
|
|
|
|
aixerr[0] += cg[0].pat[j].ixde[0];
|
|
aixerr[1] += cg[0].pat[j].ixde[1];
|
|
aixerr[2] += cg[0].pat[j].ixde[2];
|
|
|
|
if (verb >= 2) {
|
|
|
|
printf("%s%s%s: %f %f %f <=> %f %f %f de %f\n",
|
|
cg[0].pat[j].sid,
|
|
cg[0].pat[j].loc[0] != '\000' ? " " : "",
|
|
cg[0].pat[j].loc,
|
|
cg[0].pat[j].v[0], cg[0].pat[j].v[1], cg[0].pat[j].v[2],
|
|
cg[1].pat[match[j]].v[0], cg[1].pat[match[j]].v[1], cg[1].pat[match[j]].v[2],
|
|
de);
|
|
|
|
#ifdef NEVER /* Print XYZ as well */
|
|
printf(" %f %f %f <=> %f %f %f\n",
|
|
cg[0].pat[j].xyz[0], cg[0].pat[j].xyz[1], cg[0].pat[j].xyz[2],
|
|
cg[1].pat[match[j]].xyz[0], cg[1].pat[match[j]].xyz[1], cg[1].pat[match[j]].xyz[2]);
|
|
#endif
|
|
}
|
|
|
|
if (de > merr)
|
|
merr = de;
|
|
|
|
if (dovrml) {
|
|
if ((dovrml == 3 || dovrml == 4)) { /* RGB/YCC device plot */
|
|
double *val1, *val2;
|
|
int k;
|
|
|
|
if (dovrml == 3) {
|
|
val1 = cg[0].pat[i].rgb;
|
|
val2 = cg[1].pat[match[i]].rgb;
|
|
} else {
|
|
val1 = cg[0].pat[i].ycc;
|
|
val2 = cg[1].pat[match[i]].ycc;
|
|
}
|
|
|
|
de = icmNorm33(val1, val2);
|
|
|
|
if (de > 1e-6) {
|
|
wrl->add_vertex(wrl, 0, val1);
|
|
wrl->add_vertex(wrl, 0, val2);
|
|
}
|
|
|
|
#ifdef NEVER // Green target
|
|
wrl->add_marker(wrl, val1, green, rad);
|
|
|
|
#else // Natural color
|
|
for (k = 0; k < 3; k++)
|
|
col[k] = 0.3 + 0.7 * (cg[0].pat[i].rgb[k] - min[k])/(max[k] - min[k]);
|
|
wrl->add_marker(wrl, val1, col, rad);
|
|
#endif
|
|
|
|
wrl->add_marker_trans(wrl, val2, red, 0.3, rad * 0.99);
|
|
|
|
} else { /* PCS */
|
|
if (de > 1e-6) {
|
|
wrl->add_vertex(wrl, 0, cg[0].pat[i].v);
|
|
wrl->add_vertex(wrl, 0, cg[1].pat[match[i]].v);
|
|
}
|
|
if (dovrml == 2) {
|
|
wrl->add_marker(wrl, cg[0].pat[i].v, green, rad);
|
|
wrl->add_marker(wrl, cg[1].pat[match[i]].v, red, rad);
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
if (cg[0].nig > 0) {
|
|
aerr /= (double)cg[0].nig;
|
|
aierr[0] /= (double)cg[0].nig;
|
|
aierr[1] /= (double)cg[0].nig;
|
|
aierr[2] /= (double)cg[0].nig;
|
|
|
|
aixerr[0] /= (double)cg[0].nig;
|
|
aixerr[1] /= (double)cg[0].nig;
|
|
aixerr[2] /= (double)cg[0].nig;
|
|
}
|
|
|
|
if (dovrml) {
|
|
wrl->make_lines(wrl, 0, 2);
|
|
wrl->del(wrl);
|
|
wrl = NULL;
|
|
}
|
|
|
|
/* Do best 90% */
|
|
n90 = (int)(cg[0].nig * 9.0/10.0 + 0.5);
|
|
for (i = j = 0; i < cg[0].npat; i++) {
|
|
double de = cg[0].pat[sort[i]].de;
|
|
if (cg[0].pat[i].og) /* Skip out of gamut */
|
|
continue;
|
|
if (j >= (cg[0].nig-n90)) { /* If in top 90% of in gamut patches */
|
|
aerr90 += de;
|
|
if (de > merr90)
|
|
merr90 = de;
|
|
}
|
|
j++; /* Index of within gamut patches */
|
|
}
|
|
if (n90 > 0)
|
|
aerr90 /= (double)n90;
|
|
|
|
/* Do worst 10% */
|
|
n10 = (int)(cg[0].nig * 1.0/10.0 + 0.5);
|
|
for (i = j = 0; i < cg[0].npat; i++) {
|
|
double de = cg[0].pat[sort[i]].de;
|
|
if (cg[0].pat[i].og) /* Skip out of gamut */
|
|
continue;
|
|
if (j < n10) { /* If in worst 10% of in gamut patches */
|
|
aerr10 += de;
|
|
if (de > merr10)
|
|
merr10 = de;
|
|
}
|
|
j++;
|
|
}
|
|
if (n10 > 0)
|
|
aerr10 /= (double)n10;
|
|
|
|
if (verb) {
|
|
fprintf(verbo,"No of test patches in worst 10%% are = %d\n",n10);
|
|
fprintf(verbo,"No of test patches in best 90%% are = %d\n",n90);
|
|
}
|
|
|
|
/* Single number report */
|
|
if (dozrep != 0) {
|
|
if (dozrep == 1) {
|
|
printf(" %f\t", aerr);
|
|
} else if (dozrep == 2) {
|
|
printf(" %f\t", merr);
|
|
}
|
|
fflush(stdout);
|
|
|
|
} else {
|
|
printf("Verify results:\n");
|
|
if (norm == 4)
|
|
printf(" L*a*b* ref. = average XYZ %f %f %f\n",cg[0].w[0],cg[0].w[1],cg[0].w[2]);
|
|
else if (norm == 1) {
|
|
printf(" File 1 L* ref. Y %f\n", cg[0].w[1]);
|
|
printf(" File 2 L* ref. Y %f\n", cg[1].w[1]);
|
|
} else if (norm == 2) {
|
|
printf(" File 1 L*a*b* ref. XYZ %f %f %f\n", cg[0].w[0],cg[0].w[1],cg[0].w[2]);
|
|
printf(" File 2 L*a*b* ref. XYZ %f %f %f\n", cg[1].w[0],cg[1].w[1],cg[1].w[2]);
|
|
} else if (norm == 3) {
|
|
printf(" File 1 L* ref. X+Y+Z %f %f %f\n", cg[0].w[0],cg[0].w[1],cg[0].w[2]);
|
|
printf(" File 2 L* ref. X+Y+Z %f %f %f\n", cg[1].w[0],cg[1].w[1],cg[1].w[2]);
|
|
}
|
|
printf(" Total errors%s: peak = %f, avg = %f\n", cie2k ? " (CIEDE2000)" : cie94 ? " (CIE94)" : "", merr, aerr);
|
|
printf(" Worst 10%% errors%s: peak = %f, avg = %f\n", cie2k ? " (CIEDE2000)" : cie94 ? " (CIE94)" : "", merr10, aerr10);
|
|
printf(" Best 90%% errors%s: peak = %f, avg = %f\n", cie2k ? " (CIEDE2000)" : cie94 ? " (CIE94)" : "", merr90, aerr90);
|
|
printf(" avg err X %f, Y %f, Z %f\n", aixerr[0], aixerr[1], aixerr[2]);
|
|
printf(" avg err L* %f, a* %f, b* %f\n", aierr[0], aierr[1], aierr[2]);
|
|
}
|
|
|
|
free(sort);
|
|
free(match);
|
|
free(cg[0].pat);
|
|
free(cg[1].pat);
|
|
}
|
|
|
|
if (luo != NULL)
|
|
luo->del(luo);
|
|
if (xicco != NULL)
|
|
xicco->del(xicco); /* Expansion wrapper */
|
|
if (icco != NULL)
|
|
icco->del(icco); /* Icc */
|
|
if (fp != NULL)
|
|
fp->del(fp);
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
|
|
|