69
69
/* Illuminant spectra */
71
71
/* Dummy "no illuminant" illuminant spectra used to signal an emmission */
72
/* ore equal energy 'E' illuminant */
72
/* or equal energy 'E' illuminant */
73
73
static xspect il_none = {
74
74
54, 300.0, 830.0, /* 54 bands from 300 to 830 in 10nm steps */
75
75
1.0, /* Scale factor */
1978
1977
#endif /* STOCKFWA */
1981
/* save a spectrum to a CGATS file */
1979
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1980
/* save a set of spectrum to a CGATS file */
1981
/* type 0 = SPECT, 1 = CMF */
1982
1982
/* Return NZ on error */
1983
int write_xspect(char *fname, xspect *sp) {
1983
int write_nxspect(char *fname, xspect *sp, int nspec, int type) {
1985
1985
time_t clk = time(0);
1986
1986
struct tm *tsp = localtime(&clk);
1987
1987
char *atm = asctime(tsp); /* Ascii time */
1988
1988
cgats *ocg; /* output cgats structure */
1989
1989
cgats_set_elem *setel; /* Array of set value elements */
1992
1992
/* Setup output cgats file */
1993
1993
ocg = new_cgats(); /* Create a CGATS structure */
1994
ocg->add_other(ocg, "SPECT"); /* our special type is spectral power or reflectance */
1995
ocg->add_other(ocg, "CMF"); /* our special type is spectral color matching func */
1997
ocg->add_other(ocg, "SPECT"); /* our special type is spectral power or reflectance */
1995
1998
ocg->add_table(ocg, tt_other, 0); /* Start the first table */
1997
2000
ocg->add_kword(ocg, 0, "DESCRIPTOR", "Argyll Spectral power/reflectance information",NULL);
2018
2021
ocg->add_field(ocg, 0, buf, r_t);
2021
if ((setel = (cgats_set_elem *)malloc(sizeof(cgats_set_elem) * sp->spec_n)) == NULL)
2024
if ((setel = (cgats_set_elem *)malloc(sizeof(cgats_set_elem) * sp->spec_n)) == NULL) {
2024
for (i = 0; i < sp->spec_n; i++) {
2025
setel[i].d = sp->spec[i];
2028
ocg->add_setarr(ocg, 0, setel);
2029
for (j = 0; j < nspec; j++) {
2030
for (i = 0; i < sp[j].spec_n; i++) {
2031
setel[i].d = sp[j].spec[i];
2033
ocg->add_setarr(ocg, 0, setel);
2030
2036
if (ocg->write_name(ocg, fname)) {
2043
/* restore a spectrum from a CGATS file */
2049
/* restore a set of spectrum from a CGATS file. */
2050
/* Up to nspec will be restored starting at offset off.. */
2051
/* The number restored from the file will be written to *nret */
2052
/* type = any, 1 = SPECT, 2 = CMF, 3 = both */
2044
2053
/* Return NZ on error */
2045
2054
/* (Would be nice to return an error message!) */
2046
int read_xspect(xspect *sp, char *fname) {
2055
int read_nxspect(xspect *sp, char *fname, int *nret, int off, int nspec, int type) {
2047
2056
cgats *icg; /* input cgats structure */
2058
int sflds[XSPECT_MAX_BANDS];
2051
2062
/* Open and look at the spectrum file */
2052
2063
if ((icg = new_cgats()) == NULL) { /* Create a CGATS structure */
2054
2065
printf("new_cgats() failed");
2058
icg->add_other(icg, "SPECT"); /* our special input type is spectral power or reflectance */
2071
icg->add_other(icg, ""); /* Allow any signature file */
2074
icg->add_other(icg, "SPECT"); /* Allow any signature file */
2076
icg->add_other(icg, "CMF"); /* Allow any signature file */
2060
2079
if (icg->read_name(icg, fname)) {
2062
2081
printf("CGATS file read error : %s\n",icg->err);
2067
if (icg->ntables == 0 || icg->t[0].tt != tt_other || icg->t[0].oi != 0) {
2069
printf ("Input file isn't a SPECT format file\n");
2073
2087
if (icg->ntables != 1) {
2075
2089
printf ("Input file doesn't contain exactly one table\n");
2082
2097
printf ("Input file doesn't contain keyword SPECTRAL_BANDS\n");
2086
sp->spec_n = atoi(icg->t[0].kdata[ii]);
2102
proto.spec_n = atoi(icg->t[0].kdata[ii]);
2087
2103
if ((ii = icg->find_kword(icg, 0, "SPECTRAL_START_NM")) < 0) {
2089
2105
printf ("Input file doesn't contain keyword SPECTRAL_START_NM\n");
2093
sp->spec_wl_short = atof(icg->t[0].kdata[ii]);
2110
proto.spec_wl_short = atof(icg->t[0].kdata[ii]);
2094
2111
if ((ii = icg->find_kword(icg, 0, "SPECTRAL_END_NM")) < 0) {
2096
2113
printf ("Input file doesn't contain keyword SPECTRAL_END_NM\n");
2100
sp->spec_wl_long = atof(icg->t[0].kdata[ii]);
2118
proto.spec_wl_long = atof(icg->t[0].kdata[ii]);
2102
2120
if ((ii = icg->find_kword(icg, 0, "SPECTRAL_NORM")) < 0) {
2104
2122
printf ("Input file doesn't contain keyword SPECTRAL_NORM\n");
2108
sp->norm = atof(icg->t[0].kdata[ii]);
2127
proto.norm = atof(icg->t[0].kdata[ii]);
2110
2129
/* Find the fields for spectral values */
2111
for (j = 0; j < sp->spec_n; j++) {
2130
for (i = 0; i < proto.spec_n; i++) {
2114
2133
/* Compute nearest integer wavelength */
2115
nm = (int)(XSPECT_XWL(sp, j) + 0.5);
2134
nm = (int)(XSPECT_XWL(&proto, i) + 0.5);
2116
2135
sprintf(buf,"SPEC_%03d",nm);
2118
2137
if ((fi = icg->find_field(icg, 0, buf)) < 0) {
2120
2139
printf("Input file doesn't contain field %s\n",buf);
2127
2147
printf ("Field %s in specrum is wrong type - should be a float\n",buf);
2132
sp->spec[j] = *((double *)icg->t[0].fdata[0][fi]);
2155
/* Read all the spectra */
2156
for (j = off; j < nspec && j < icg->t[0].nsets; j++) {
2158
XSPECT_COPY_INFO(&sp[j], &proto);
2160
for (i = 0; i < proto.spec_n; i++) {
2161
sp[j].spec[i] = *((double *)icg->t[0].fdata[j][sflds[i]]);
2172
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2174
/* save a spectrum to a CGATS file */
2175
/* Return NZ on error */
2176
int write_xspect(char *fname, xspect *sp) {
2177
return write_nxspect(fname, sp, 1, 0);
2180
/* restore a spectrum from a CGATS file */
2181
/* Return NZ on error */
2182
/* (Would be nice to return an error message!) */
2183
int read_xspect(xspect *sp, char *fname) {
2186
if ((rv = read_nxspect(sp, fname, &nret, 0, 1, 1)) != 0)
2190
printf ("Didn't read one spectra\n");
2198
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2199
/* save a set of 3 spectrum to a CGATS CMF file */
2200
/* Return NZ on error */
2201
int write_cmf(char *fname, xspect sp[3]) {
2202
return write_nxspect(fname, sp, 3, 1);
2205
/* restore a spectrum from a CGATS file */
2206
/* Return NZ on error */
2207
/* (Would be nice to return an error message!) */
2208
int read_cmf(xspect sp[3], char *fname) {
2211
if ((rv = read_nxspect(sp, fname, &nret, 0, 3, 2)) != 0)
2215
printf ("Didn't read three spectra\n");
2139
2224
#endif /* !SALONEINSTLIB */
2141
2227
/* Get a raw 3rd order polinomial interpolated spectrum value. */
2142
2228
/* Return NZ if value is valid, Z and last valid value */
2143
2229
/* if outside the range */
2953
3039
/* Do the FWA corrected spectral to CIE conversion. */
2954
3040
/* Note that the input spectrum normalisation value is used. */
3041
/* Emissive spectral values are assumed to be in mW/nm, and sampled */
3042
/* rather than integrated if they are not at 1nm spacing. */
2955
3043
static void xsp2cie_fwa_sconvert(
2956
3044
xsp2cie *p, /* this */
2957
3045
xspect *sout, /* Return corrected input spectrum (may be NULL, or same as imput) */
2982
3070
tsout.spec_wl_long = 0.0;
2983
3071
tsout.norm = 0.0;
3073
#define MIN_UVILLUM 1e-8 /* Minimum assumed UV illumination level at wavelength */
3074
#define MIN_UVREFL 1e-6 /* Minimum assumed UV reflectance at wavelength */
3075
#define MIN_ILLUM 1e-6 /* Minimum assumed illumination level at wavelength */
3076
#define MIN_REFL 1e-6 /* Minimum assumed reflectance at wavelength */
2985
3078
/* With colorant, estimate stimulation level of FWA for instrument illuminant */
2986
3079
/* and for target illuminant. Because the colorant estimate depends on the FWA */
2987
3080
/* estimate, and the FWA emissions can contribute to FWA stimulation, */
3005
3098
Kct = Emct * Eu; /* FWA contribution under target illum. */
3007
3100
getval_lxspec(&p->instr, &Ii, ww); /* Normalised instr. illuminant at wavelength */
3101
if (Ii < MIN_UVILLUM)
3010
3103
getval_lxspec(&p->illum, &It, ww); /* Normalised target. illuminant at wavelength */
3104
if (It < MIN_UVILLUM)
3013
3106
getval_lxspec(in, &Rc, ww) ; /* Media + colorant reflectance at wavelength */
3015
3108
getval_lxspec(&p->media, &Rmb, ww); /* Base media reflectance at this wavelength */
3109
if (Rmb < MIN_UVREFL)
3020
3113
Rcch = sqrt(Rc/Rmb); /* Half reflectance estimate (valid if no FWA) */
3023
3116
/* Solve for underlying colorant half reflectance, discounting FWA */
3024
if (Rmb < 1e-9) /* Hmm. */
3117
if (Rmb <= MIN_UVREFL) /* Hmm. */
3025
3118
Rcch = sqrt(fabs(Rmb));
3027
3120
Rcch = (-Kc + sqrt(Kc * Kc + 4.0 * Ii * Ii * Rmb * Rc))/(2.0 * Ii * Rmb);
3030
3123
getval_lxspec(&FWA1_stim, &Su, ww); /* FWA stimulation sensitivity this wavelength */
3031
3126
Smc += Su * (Ii * Rcch + Kc);
3032
3127
Smct += Su * (It * Rcch + Kct);
3129
printf("~1 at %.1fnm, Rmb %f, Rc %f, Rch %f, Rcch %f, Ii %f, It %f, Kct %f, Smc %f, Smct %f,\n",ww,Rmb,Rc,sqrt(Rc),Rcch,Ii,It,Kct,Su * (Ii * Rcch + Kc),Su * (It * Rcch + Kct));
3034
3132
Emc = Smc/p->Sm; /* FWA Emmsion muliplier with colorant for instr. illum. */
3035
3133
Emct = Smct/p->Sm; /* FWA Emmsion muliplier with colorant for target illum. */
3038
printf("~1 Smc = %f\n",Smc); fflush(stdout);
3039
printf("~1 Smct = %f\n",Smct); fflush(stdout);
3040
printf("~1 Emc = %f\n",Emc); fflush(stdout);
3041
printf("~1 Emct = %f\n",Emct); fflush(stdout);
3136
printf("~1 Itteration %d, Smc %f, Smct %f, Emc %f, Emct %f\n",k, Smc,Smct,Emc,Emct); fflush(stdout);
3059
3154
double Rmb; /* Base media reflectance estimate */
3060
3155
double Eu; /* FWA emmission profile */
3061
3156
double Rc; /* Measured reflectance under inst. illum. */
3157
/* Rch Measured half reflectance under inst. illum */
3062
3158
double Rcch; /* Corrected Rc colorant half reflectance */
3063
3159
double RctI; /* Corrected Rc for target illuminant times illuminant */
3069
3165
getval_lxspec(&p->media, &Rmb, ww); /* Base Media */
3070
3166
getval_lxspec(in, &Rc, ww); /* Media + colorant reflectance at wavelength + FWA */
3071
3167
getval_lxspec(&p->instr, &Ii, ww); /* Normalised instrument illuminant */
3075
/* Solve for underlying colorant half reflectance, discounting FWA */
3076
if (Rmb < 1e-9) /* Hmm. */
3171
/* Solve for underlying colorant half transmittance, discounting FWA */
3172
if (Rmb <= MIN_REFL) /* Hmm. */
3077
3173
Rcch = sqrt(fabs(Rmb));
3079
3175
Rcch = (-Kc + sqrt(Kc * Kc + 4.0 * Ii * Ii * Rmb * Rc))/(2.0 * Ii * Rmb);
3082
printf("~1 at %fnm, Rc = %f, Rch = %f, Rcch = %f\n",ww,Rc,sqrt(Rc),Rcch);
3084
3177
/* Estimated reflectance times target illum. */
3085
3178
getval_lxspec(&p->illum, &It, ww); /* Normalised target illuminant */
3088
3181
RctI = (It * Rcch * Rmb + Kct) * Rcch;
3184
printf("~1 at %.1fnm, Rmb %f, Rc %f, Rch %f, Rcch %f, Ii %f, It %f, Kct %f, RctI %f, CrdRef %f\n",ww,Rmb,Rc,sqrt(Rc),Rcch,Ii,It,Kct,RctI,RctI/It);
3090
3192
#ifdef DOPLOT_ALL_FWA
3092
3194
y1[plix] = Rc; /* Uncorrected reflectance */
3433
3536
/* Do the normal spectral to CIE conversion. */
3434
3537
/* Note that the input spectrum normalisation value is used. */
3538
/* Emissive spectral values are assumed to be in mW/nm, and sampled */
3539
/* rather than integrated if they are not at 1nm spacing. */
3435
3540
void xsp2cie_sconvert(
3436
3541
xsp2cie *p, /* this */
3437
3542
xspect *sout, /* Return input spectrum (may be NULL) */
3448
3553
/* Integrate at 1nm intervals over the observer range (as */
3449
3554
/* per CIE recommendations). Lower resolution spectra are */
3450
/* upsampled using linear/3rd order poolinomial interpolated */
3555
/* upsampled using linear/3rd order polinomial interpolated */
3451
3556
/* (also as per CIE recommendations), and consistent (?) with the */
3452
3557
/* assumption of a triangular spectral response made in the */
3453
3558
/* ANSI CGATS.5-1993 spec. If illumninant or material spectra */
3490
/* Normal conversion */
3596
/* Normal Tristumulus conversion */
3491
3597
void xsp2cie_convert(xsp2cie *p, double *out, xspect *in) {
3492
3598
xsp2cie_sconvert(p, NULL, out, in);
3504
3610
icxIllumeType ilType, /* Illuminant */
3505
3611
xspect *custIllum, /* Optional custom illuminant */
3506
3612
icxObserverType obType, /* Observer */
3507
xspect *custObserver[3] /* Optional custom observer */
3508
#ifndef SALONEINSTLIB
3509
, icColorSpaceSignature rcs /* Return color space, icSigXYZData or icSigLabData */
3510
#endif /* !SALONEINSTLIB */
3613
xspect custObserver[3], /* Optional custom observer */
3614
icColorSpaceSignature rcs /* Return color space, icSigXYZData or icSigLabData */
3615
/* ** Must be icSigXYZData if SALONEINSTLIB ** */
3564
3669
/* Do 3 structure copies to record observer sensitivity curves */
3565
3670
switch (obType) {
3566
3671
case icxOT_custom:
3567
p->observer[0] = *custObserver[0];
3568
p->observer[1] = *custObserver[1];
3569
p->observer[2] = *custObserver[2];
3672
p->observer[0] = custObserver[0];
3673
p->observer[1] = custObserver[1];
3674
p->observer[2] = custObserver[2];
3571
3676
case icxOT_default:
3572
3677
case icxOT_CIE_1931_2:
4017
4122
int icx_ill_sp2XYZ(
4018
4123
double xyz[3], /* Return XYZ value with Y == 1 */
4019
4124
icxObserverType obType, /* Observer */
4020
xspect *custObserver[3],/* Optional custom observer */
4125
xspect custObserver[3], /* Optional custom observer */
4021
4126
icxIllumeType ilType, /* Type of illuminant, icxIT_Dtemp or icxIT_Ptemp */
4022
4127
double ct, /* Input temperature in degrees K */
4023
4128
xspect *custIllum /* Optional custom illuminant */
4107
4212
double txyz[3], /* If not NULL, return the XYZ of the locus temperature */
4108
4213
icxIllumeType ilType, /* Type of illuminant, icxIT_Dtemp or icxIT_Ptemp */
4109
4214
icxObserverType obType, /* Observer */
4110
xspect *custObserver[3],/* Optional custom observer */
4215
xspect custObserver[3], /* Optional custom observer */
4111
4216
double xyz[3], /* Input XYZ value, NULL if spectrum intead */
4112
4217
xspect *insp, /* Input spectrum value, NULL if xyz[] instead */
4113
4218
int viscct /* nz to use visual CIEDE2000, 0 to use CCT CIE 1960 UCS. */
added
removed
xicc/xspect.c
