541 lines
13 KiB
C
541 lines
13 KiB
C
|
|
/*
|
|
* International Color Consortium color transform expanded support
|
|
*
|
|
* Author: Graeme W. Gill
|
|
* Date: 2006/5/9
|
|
* Version: 1.00
|
|
*
|
|
* Copyright 2006 Graeme W. Gill
|
|
* All rights reserved.
|
|
* This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
|
|
* see the License.txt file for licencing details.
|
|
*
|
|
*/
|
|
|
|
/*
|
|
* This is some test code to test the Daylight and Plankian spectra,
|
|
* Correlated and Visual Color Temperatures, and CRI.
|
|
* and plot a spectrum, CMF or CCSS.
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <stdarg.h>
|
|
#include <math.h>
|
|
#include "aconfig.h"
|
|
#include "cgats.h"
|
|
#include "xspect.h"
|
|
#include "numlib.h"
|
|
#include "plot.h"
|
|
#include "ui.h"
|
|
|
|
#define PLANKIAN
|
|
#define XRES 500
|
|
|
|
#define MAXGRAPHS 10
|
|
|
|
/* Display a spectrum etc. */
|
|
/* We are guaranteed that the x range/increments are identical, */
|
|
/* and that there is only one spectrum if douv */
|
|
static int do_spec(
|
|
char name[MAXGRAPHS][200],
|
|
xspect *sp,
|
|
int nsp, /* Number of sp */
|
|
inst_meas_type mt, /* Measurement type */
|
|
int dozero, /* Include zero in the range */
|
|
int noplot, /* Don't actually plot */
|
|
icxDensityType dens, /* Density type to print */
|
|
int douv, /* Do variation of added UV test */
|
|
icxObserverType obType, /* Observer (defaults to 1931 2 degree) */
|
|
xspect *custObserver, /* If obType = icxOT_custom */
|
|
double uvmin,
|
|
double uvmax
|
|
) {
|
|
int n, i, j, k, m;
|
|
double wl_short, wl_long; /* Common range */
|
|
double xyz[3]; /* Color temperature */
|
|
double Yxy[3];
|
|
double Lab[3]; /* D50 Lab value */
|
|
double xx[XRES];
|
|
double yy[10][XRES];
|
|
double *yp[10];
|
|
double cct, vct;
|
|
double cct_xyz[3], vct_xyz[3];
|
|
double cct_lab[3], vct_lab[3];
|
|
double cdt, vdt;
|
|
double cdt_xyz[3], vdt_xyz[3];
|
|
double cdt_lab[3], vdt_lab[3];
|
|
icmXYZNumber wp;
|
|
double de;
|
|
double uv = uvmin;
|
|
double step = 0.1;
|
|
xspect tsp; /* Spectrum with possible UV added */
|
|
char *color[] = {
|
|
"Black", "Red", "Green", "Blue", "Yellow", "Purple", "Brown", "Orange", "Grey", "Magenta"
|
|
};
|
|
|
|
printf("\n");
|
|
|
|
for (j = 0; j < 10; j++)
|
|
yp[j] = NULL;
|
|
|
|
if (nsp > 10)
|
|
nsp = 10;
|
|
|
|
m = 0; /* offset in output array */
|
|
n = 1;
|
|
|
|
wl_short = 1e6;
|
|
wl_long = -1e6;
|
|
for (k = 0; k < nsp; k++) {
|
|
if (sp[k].spec_wl_long > wl_long)
|
|
wl_long = sp[k].spec_wl_long;
|
|
if (sp[k].spec_wl_short < wl_short)
|
|
wl_short = sp[k].spec_wl_short;
|
|
}
|
|
|
|
if (douv) {
|
|
n = 1 + (int)(0.5 + (uvmax-uvmin)/0.1);
|
|
if (n > 9)
|
|
n = 9; /* Don't use white */
|
|
if (n > 1)
|
|
step = (uvmax-uvmin)/(n-1.0);
|
|
}
|
|
|
|
for (k = 0; k < nsp; k++) {
|
|
tsp = sp[k];
|
|
for (uv = uvmax, j = 0; j < n; j++, uv -= step) {
|
|
|
|
if (douv) {
|
|
printf("UV level = %f\n",uv);
|
|
xsp_setUV(&tsp, &sp[k], uv);
|
|
}
|
|
|
|
printf("Type = %s [%s]\n",name[k], color[k]);
|
|
|
|
if (mt == inst_mrt_none
|
|
|| mt == inst_mrt_emission
|
|
|| mt == inst_mrt_ambient
|
|
|| mt == inst_mrt_emission_flash
|
|
|| mt == inst_mrt_ambient_flash) {
|
|
|
|
/* Compute absolute XYZ of illuminant */
|
|
if (icx_ill_sp2XYZ(xyz, obType, custObserver, icxIT_custom, 0, &tsp, 1) != 0)
|
|
warning("icx_ill_sp2XYZ returned error");
|
|
|
|
printf("Abs. Y = %f\n", xyz[1]);
|
|
|
|
/* Normalise to Y = 1 */
|
|
xyz[0] /= xyz[1];
|
|
xyz[2] /= xyz[1];
|
|
xyz[1] /= xyz[1];
|
|
|
|
icmXYZ2Yxy(Yxy, xyz);
|
|
icmXYZ2Lab(&icmD50, Lab, xyz);
|
|
|
|
printf("XYZ = %f %f %f, x,y = %f %f\n", xyz[0], xyz[1], xyz[2], Yxy[1], Yxy[2]);
|
|
printf("D50 L*a*b* = %f %f %f\n", Lab[0], Lab[1], Lab[2]);
|
|
|
|
/* Compute CCT */
|
|
if ((cct = icx_XYZ2ill_ct(cct_xyz, icxIT_Ptemp, obType, custObserver, xyz, NULL, 0)) < 0)
|
|
warning("Got bad cct\n");
|
|
|
|
/* Compute VCT */
|
|
if ((vct = icx_XYZ2ill_ct(vct_xyz, icxIT_Ptemp, obType, custObserver, xyz, NULL, 1)) < 0)
|
|
warning("Got bad vct\n");
|
|
|
|
printf("CCT = %f, VCT = %f\n",cct, vct);
|
|
|
|
/* Compute CDT */
|
|
if ((cct = icx_XYZ2ill_ct(cct_xyz, icxIT_Dtemp, obType, custObserver, xyz, NULL, 0)) < 0)
|
|
warning("Got bad cct\n");
|
|
|
|
/* Compute VDT */
|
|
if ((vct = icx_XYZ2ill_ct(vct_xyz, icxIT_Dtemp, obType, custObserver, xyz, NULL, 1)) < 0)
|
|
warning("Got bad vct\n");
|
|
|
|
printf("CDT = %f, VDT = %f\n",cct, vct);
|
|
|
|
{
|
|
int invalid = 0;
|
|
double RR[14];
|
|
double cri;
|
|
cri = icx_CIE1995_CRI(&invalid, RR, &tsp);
|
|
printf("CRI = %.1f [ R9 = %.1f ]%s\n",cri,RR[9-1],invalid ? " (Invalid)" : "");
|
|
}
|
|
{
|
|
int invalid = 0;
|
|
double tlci;
|
|
tlci = icx_EBU2012_TLCI(&invalid, &tsp);
|
|
printf("TLCI = %.1f%s\n",tlci,invalid ? " (Invalid)" : "");
|
|
}
|
|
|
|
/* Use modern color difference - gives a better visual match */
|
|
icmAry2XYZ(wp, vct_xyz);
|
|
icmXYZ2Lab(&wp, cct_lab, cct_xyz);
|
|
icmXYZ2Lab(&wp, vct_lab, vct_xyz);
|
|
de = icmCIE2K(cct_lab, vct_lab);
|
|
printf("CIEDE2000 Delta E = %f\n",de);
|
|
|
|
} else if (mt == inst_mrt_none
|
|
|| mt == inst_mrt_reflective
|
|
|| mt == inst_mrt_transmissive) {
|
|
double cmyv[4];
|
|
|
|
printf("CIE values under D50 illuminant:\n");
|
|
|
|
if (icx_sp2XYZ(xyz, obType, custObserver, icxIT_D50, 0, NULL, &tsp) != 0)
|
|
warning("icx_sp2XYZ returned error");
|
|
|
|
icmXYZ2Yxy(Yxy, xyz);
|
|
icmXYZ2Lab(&icmD50, Lab, xyz);
|
|
|
|
printf("XYZ = %f %f %f, x,y = %f %f\n", xyz[0], xyz[1], xyz[2], Yxy[1], Yxy[2]);
|
|
printf("D50 L*a*b* = %f %f %f\n", Lab[0], Lab[1], Lab[2]);
|
|
|
|
if (dens != icxDT_none) {
|
|
/* Density */
|
|
xsp_density(cmyv, &tsp, dens);
|
|
printf("CMYV %s = %f %f %f %f\n", xsp_density_desc(dens), cmyv[0], cmyv[1], cmyv[2], cmyv[3]);
|
|
}
|
|
|
|
} else {
|
|
printf("Unhandled measurement type '%s'\n",meas_type2str(mt));
|
|
}
|
|
|
|
/* Plot spectrum out */
|
|
for (i = 0; i < XRES; i++) {
|
|
double ww;
|
|
|
|
ww = (wl_long - wl_short)
|
|
* ((double)i/(XRES-1.0)) + wl_short;
|
|
|
|
xx[i] = ww;
|
|
yy[(m + k + j) % 10][i] = value_xspect(&tsp, ww);
|
|
}
|
|
yp[(m + k + j) % 10] = &yy[(m + k + j) % 10][0];
|
|
printf("\n");
|
|
}
|
|
}
|
|
if (!noplot) {
|
|
do_plot10(xx, yp[0], yp[1], yp[2], yp[3], yp[4], yp[5], yp[6], yp[7], yp[8], yp[9], XRES, dozero);
|
|
}
|
|
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
void usage(void) {
|
|
fprintf(stderr,"Plot spectrum and calculate CCT and VCT\n");
|
|
fprintf(stderr,"Author: Graeme W. Gill, licensed under the AGPL Version 3\n");
|
|
fprintf(stderr,"usage: specplot [infile.sp]\n");
|
|
fprintf(stderr," -v verbose\n");
|
|
fprintf(stderr," -c combine multiple files into one plot\n");
|
|
fprintf(stderr," -n normalize max level\n");
|
|
fprintf(stderr," -z don't make range cover zero\n");
|
|
fprintf(stderr," -s don't plot spectra\n");
|
|
fprintf(stderr," -d I|A|M|T|E print density values\n");
|
|
fprintf(stderr," -u level plot effect of adding estimated UV level\n");
|
|
fprintf(stderr," -U plot effect of adding range of estimated UV level\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," [infile.sp ...] spectrum files to plot\n");
|
|
fprintf(stderr," default is all built in illuminants\n");
|
|
exit(1);
|
|
}
|
|
|
|
int main(int argc, char *argv[]) {
|
|
int fa, nfa; /* argument we're looking at */
|
|
int k;
|
|
int verb = 0;
|
|
int comb = 0;
|
|
int norm = 0;
|
|
int zero = 1;
|
|
int noplot = 0;
|
|
icxDensityType dens = icxDT_none;
|
|
double temp;
|
|
icxObserverType obType = icxOT_CIE_1931_2; /* Observer */
|
|
xspect custObserver[3]; /* If obType = icxOT_custom */
|
|
xspect sp[MAXGRAPHS];
|
|
icxIllumeType ilType;
|
|
int douv = 0;
|
|
double uvmin = -1.0, uvmax = 1.0;
|
|
char buf[MAXGRAPHS][200];
|
|
|
|
|
|
error_program = argv[0];
|
|
|
|
/* Process the arguments */
|
|
for (fa = 1;fa < argc;fa++) {
|
|
nfa = fa; /* skip to nfa if next argument is used */
|
|
if (argv[fa][0] == '-') { /* Look for any flags */
|
|
char *na = NULL; /* next argument after flag, null if none */
|
|
|
|
if (argv[fa][2] != '\000')
|
|
na = &argv[fa][2]; /* next is directly after flag */
|
|
else {
|
|
if ((fa+1) < argc) {
|
|
if (argv[fa+1][0] != '-') {
|
|
nfa = fa + 1;
|
|
na = argv[nfa]; /* next is seperate non-flag argument */
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Show added UV */
|
|
if (argv[fa][1] == 'u') {
|
|
douv = 1;
|
|
|
|
fa = nfa;
|
|
if (na == NULL)
|
|
usage();
|
|
|
|
uvmin = uvmax = atof(na);
|
|
if (uvmin < -10.0 || uvmax > 10.0)
|
|
usage();
|
|
}
|
|
|
|
else if (argv[fa][1] == 'U') {
|
|
douv = 1;
|
|
}
|
|
|
|
/* Spectral Observer type */
|
|
else if (argv[fa][1] == 'o') {
|
|
if (na == NULL) usage();
|
|
fa = nfa;
|
|
if (strcmp(na, "1931_2") == 0) { /* Classic 2 degree */
|
|
obType = icxOT_CIE_1931_2;
|
|
} else if (strcmp(na, "1964_10") == 0) { /* Classic 10 degree */
|
|
obType = icxOT_CIE_1964_10;
|
|
} else if (strcmp(na, "2015_2") == 0) { /* Latest 2 degree */
|
|
obType = icxOT_CIE_2015_2;
|
|
} else if (strcmp(na, "2015_10") == 0) { /* Latest 10 degree */
|
|
obType = icxOT_CIE_2015_10;
|
|
} else if (strcmp(na, "1955_2") == 0) { /* Stiles and Burch 1955 2 degree */
|
|
obType = icxOT_Stiles_Burch_2;
|
|
} else if (strcmp(na, "1978_2") == 0) { /* Judd and Voss 1978 2 degree */
|
|
obType = icxOT_Judd_Voss_2;
|
|
} else if (strcmp(na, "shaw") == 0) { /* Shaw and Fairchilds 1997 2 degree */
|
|
obType = icxOT_Shaw_Fairchild_2;
|
|
} else { /* Assume it's a filename */
|
|
obType = icxOT_custom;
|
|
if (read_cmf(custObserver, na) != 0) {
|
|
fprintf(stderr,"Failed to read custom observer CMF from -o file '%s'\n",na);
|
|
usage();
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Verbosity */
|
|
else if (argv[fa][1] == 'v') {
|
|
verb = 1;
|
|
|
|
} else if (argv[fa][1] == 'c') {
|
|
comb = 1;
|
|
|
|
} else if (argv[fa][1] == 'n') {
|
|
norm = 1;
|
|
|
|
} else if (argv[fa][1] == 'z') {
|
|
zero = 0;
|
|
|
|
} else if (argv[fa][1] == 's') {
|
|
noplot = 1;
|
|
|
|
} else if (argv[fa][1] == 'd') {
|
|
if (na == NULL)
|
|
usage();
|
|
|
|
switch(na[0]) {
|
|
case 'N':
|
|
dens = icxDT_none;
|
|
break;
|
|
case 'I':
|
|
dens = icxDT_ISO;
|
|
break;
|
|
case 'A':
|
|
dens = icxDT_A;
|
|
break;
|
|
case 'M':
|
|
dens = icxDT_M;
|
|
break;
|
|
case 'T':
|
|
dens = icxDT_T;
|
|
break;
|
|
case 'E':
|
|
dens = icxDT_E;
|
|
break;
|
|
default:
|
|
usage();
|
|
}
|
|
fa = nfa;
|
|
|
|
} else {
|
|
usage();
|
|
}
|
|
}
|
|
else
|
|
break;
|
|
}
|
|
|
|
|
|
if (fa < argc && argv[fa][0] != '-') { /* Got file arguments */
|
|
int nsp = 0; /* Current number in sp[] */
|
|
int soff = 0; /* Offset within file */
|
|
int maxgraphs = MAXGRAPHS;
|
|
int eof;
|
|
|
|
if (douv)
|
|
maxgraphs = 1;
|
|
|
|
nsp = 0;
|
|
|
|
/* Until we run out */
|
|
for (;;) {
|
|
int i, nret, nreq;
|
|
inst_meas_type mt;
|
|
inst_meas_cond mc;
|
|
|
|
/* If we've got to the limit of each plot, */
|
|
/* or at least one and we're not combining files and at start of a new file, */
|
|
/* or at least one and there are no more files */
|
|
if (nsp >= MAXGRAPHS || (nsp > 0 && ((!comb && soff == 0) || fa >= argc))) {
|
|
if (norm) {
|
|
int j;
|
|
double mlev = 0.0;
|
|
|
|
for (i = 0; i < nsp; i++) {
|
|
|
|
xspect_denorm(&sp[i]);
|
|
|
|
for (j = 0; j < sp[i].spec_n; j++) {
|
|
if (sp[i].spec[j] > mlev)
|
|
mlev = sp[i].spec[j];
|
|
}
|
|
}
|
|
|
|
if (mlev > 1e-6) {
|
|
for (i = 0; i < nsp; i++) {
|
|
double mm = 0.0;
|
|
|
|
for (j = 0; j < sp[i].spec_n; j++) {
|
|
if (sp[i].spec[j] > mm)
|
|
mm = sp[i].spec[j];
|
|
}
|
|
|
|
if (mm > 1e-6) {
|
|
mm = mlev/mm;
|
|
for (j = 0; j < sp[i].spec_n; j++) {
|
|
sp[i].spec[j] *= mm;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Plot what we've got */
|
|
do_spec(buf, sp, nsp, mt, zero, noplot, dens, douv, obType, custObserver, uvmin, uvmax);
|
|
nsp = 0;
|
|
}
|
|
|
|
if (fa >= argc) /* No more files */
|
|
break;
|
|
|
|
/* Read as many spectra from the file as possible */
|
|
nreq = MAXGRAPHS - nsp;
|
|
|
|
if (read_nxspect(&sp[nsp], &mt, &mc, argv[fa], &nret, soff, nreq, 0) != 0) {
|
|
error ("Unable to read custom spectrum, CMF or CCSS '%s'",argv[fa]);
|
|
}
|
|
|
|
for (i = 0; i < nret; i++) {
|
|
xspect_denorm(&sp[nsp + i]);
|
|
sprintf(buf[nsp + i],"File '%s' spect %d",argv[fa], soff + i);
|
|
}
|
|
nsp += nret;
|
|
soff += nret;
|
|
if (nret < nreq) { /* We're done with this file */
|
|
fa++;
|
|
soff = 0;
|
|
}
|
|
}
|
|
|
|
} else {
|
|
|
|
/* For each standard illuminant */
|
|
for (ilType = icxIT_A; ilType <= icxIT_F10; ilType++) {
|
|
char *inm = NULL;
|
|
|
|
switch (ilType) {
|
|
case icxIT_A:
|
|
inm = "A"; break;
|
|
case icxIT_C:
|
|
inm = "C"; break;
|
|
case icxIT_D50:
|
|
inm = "D50"; break;
|
|
case icxIT_D50M2:
|
|
inm = "D50M2"; break;
|
|
case icxIT_D55:
|
|
inm = "D55"; break;
|
|
case icxIT_D65:
|
|
inm = "D65"; break;
|
|
case icxIT_D75:
|
|
inm = "D75"; break;
|
|
case icxIT_E:
|
|
inm = "E"; break;
|
|
case icxIT_F5:
|
|
inm = "F5"; break;
|
|
case icxIT_F8:
|
|
inm = "F8"; break;
|
|
case icxIT_F10:
|
|
inm = "F10"; break;
|
|
default:
|
|
inm = "Unknown"; break;
|
|
break;
|
|
}
|
|
|
|
if (standardIlluminant(&sp[0], ilType, 0) != 0)
|
|
error ("standardIlluminant returned error for %d (%s)",ilType,inm);
|
|
|
|
strcpy(buf[0],inm);
|
|
do_spec(buf, sp, 1, inst_mrt_ambient, zero, noplot, dens, douv, obType, custObserver, uvmin, uvmax);
|
|
}
|
|
|
|
/* For each material and illuminant */
|
|
for (k = 0; k < 2; k++) {
|
|
double start, end;
|
|
|
|
if (k == 0)
|
|
start = 3400, end = 9000;
|
|
else
|
|
start = 500, end = 9000;
|
|
|
|
for (temp = start; temp <= end; temp += 500) {
|
|
|
|
ilType = k == 0 ? icxIT_Dtemp : icxIT_Ptemp;
|
|
|
|
if (standardIlluminant(&sp[0], ilType, temp) != 0)
|
|
error ("standardIlluminant returned error for %s temp %f",
|
|
k == 0 ? "Daylight" : "Plankian", temp);
|
|
|
|
sprintf(buf[0], "%s at %f", k == 0 ? "Daylight" : "Black body", temp);
|
|
|
|
do_spec(buf, sp, 1, inst_mrt_ambient, zero, noplot, dens, douv, obType, custObserver, uvmin, uvmax);
|
|
}
|
|
}
|
|
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|