ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/ccolor.c
Revision: 3.13
Committed: Thu Jan 4 16:14:50 2024 UTC (3 months, 4 weeks ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 3.12: +8 -8 lines
Log Message:
fix(mgf2rad): Made color equivalence check more stringent

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: ccolor.c,v 3.12 2024/01/04 01:55:42 greg Exp $";
3 #endif
4 /*
5 * Spectral color handling routines
6 */
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include "ccolor.h"
12
13 #undef frand
14 #define frand() (rand()*(1./(RAND_MAX+.5)))
15
16 /* Sharp primary matrix */
17 float XYZtoSharp[3][3] = {
18 { 1.2694, -0.0988, -0.1706},
19 {-0.8364, 1.8006, 0.0357},
20 { 0.0297, -0.0315, 1.0018}
21 };
22 /* inverse Sharp primary matrix */
23 float XYZfromSharp[3][3] = {
24 { 0.8156, 0.0472, 0.1372},
25 { 0.3791, 0.5769, 0.0440},
26 {-0.0123, 0.0167, 0.9955}
27 };
28
29 const C_COLOR c_dfcolor = C_DEFCOLOR;
30
31 const C_CHROMA c_dfchroma = 49750; /* c_encodeChroma(&c_dfcolor) */
32
33 /* CIE 1931 Standard Observer curves */
34 const C_COLOR c_x31 = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
35 {14,42,143,435,1344,2839,3483,3362,2908,1954,956,
36 320,49,93,633,1655,2904,4334,5945,7621,9163,10263,
37 10622,10026,8544,6424,4479,2835,1649,874,468,227,
38 114,58,29,14,7,3,2,1,0}, 106836L, .467, .368, 362.230
39 };
40 const C_COLOR c_y31 = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
41 {0,1,4,12,40,116,230,380,600,910,1390,2080,3230,
42 5030,7100,8620,9540,9950,9950,9520,8700,7570,6310,
43 5030,3810,2650,1750,1070,610,320,170,82,41,21,10,
44 5,2,1,1,0,0}, 106856L, .398, .542, 493.525
45 };
46 const C_COLOR c_z31 = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
47 {65,201,679,2074,6456,13856,17471,17721,16692,
48 12876,8130,4652,2720,1582,782,422,203,87,39,21,17,
49 11,8,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
50 106770L, .147, .077, 54.363
51 };
52 /* Derived CIE 1931 Primaries (imaginary) */
53 static const C_COLOR cie_xp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
54 {-174,-198,-195,-197,-202,-213,-235,-272,-333,
55 -444,-688,-1232,-2393,-4497,-6876,-6758,-5256,
56 -3100,-815,1320,3200,4782,5998,6861,7408,7754,
57 7980,8120,8199,8240,8271,8292,8309,8283,8469,
58 8336,8336,8336,8336,8336,8336},
59 127424L, 1., .0,
60 };
61 static const C_COLOR cie_yp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
62 {-451,-431,-431,-430,-427,-417,-399,-366,-312,
63 -204,57,691,2142,4990,8810,9871,9122,7321,5145,
64 3023,1123,-473,-1704,-2572,-3127,-3474,-3704,
65 -3846,-3927,-3968,-3999,-4021,-4038,-4012,-4201,
66 -4066,-4066,-4066,-4066,-4066,-4066},
67 -23035L, .0, 1.,
68 };
69 static const C_COLOR cie_zp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
70 {4051,4054,4052,4053,4054,4056,4059,4064,4071,
71 4074,4056,3967,3677,2933,1492,313,-440,-795,
72 -904,-918,-898,-884,-869,-863,-855,-855,-851,
73 -848,-847,-846,-846,-846,-845,-846,-843,-845,
74 -845,-845,-845,-845,-845},
75 36057L, .0, .0,
76 };
77
78
79 /* convert to sharpened RGB color for low-error operations */
80 void
81 c_toSharpRGB(C_COLOR *cin, double cieY, float cout[3])
82 {
83 double xyz[3];
84
85 c_ccvt(cin, C_CSXY);
86
87 xyz[0] = cin->cx/cin->cy * cieY;
88 xyz[1] = cieY;
89 xyz[2] = (1. - cin->cx - cin->cy)/cin->cy * cieY;
90
91 cout[0] = XYZtoSharp[0][0]*xyz[0] + XYZtoSharp[0][1]*xyz[1] +
92 XYZtoSharp[0][2]*xyz[2];
93 cout[1] = XYZtoSharp[1][0]*xyz[0] + XYZtoSharp[1][1]*xyz[1] +
94 XYZtoSharp[1][2]*xyz[2];
95 cout[2] = XYZtoSharp[2][0]*xyz[0] + XYZtoSharp[2][1]*xyz[1] +
96 XYZtoSharp[2][2]*xyz[2];
97 }
98
99 /* convert back from sharpened RGB color */
100 double
101 c_fromSharpRGB(float cin[3], C_COLOR *cout)
102 {
103 double xyz[3], sf;
104
105 xyz[1] = XYZfromSharp[1][0]*cin[0] + XYZfromSharp[1][1]*cin[1] +
106 XYZfromSharp[1][2]*cin[2];
107 if (xyz[1] <= 1e-6) {
108 *cout = c_dfcolor; /* punting, here... */
109 return xyz[1];
110 }
111 xyz[0] = XYZfromSharp[0][0]*cin[0] + XYZfromSharp[0][1]*cin[1] +
112 XYZfromSharp[0][2]*cin[2];
113 xyz[2] = XYZfromSharp[2][0]*cin[0] + XYZfromSharp[2][1]*cin[1] +
114 XYZfromSharp[2][2]*cin[2];
115
116 sf = 1./(xyz[0] + xyz[1] + xyz[2]);
117
118 cout->cx = xyz[0] * sf;
119 cout->cy = xyz[1] * sf;
120 cout->flags = C_CDXY|C_CSXY;
121
122 return(xyz[1]);
123 }
124
125 /* assign arbitrary spectrum and return Y value */
126 double
127 c_sset(C_COLOR *clr, double wlmin, double wlmax, const float spec[], int nwl)
128 {
129 double yval, scale;
130 float va[C_CNSS];
131 int i, pos, n, imax, wl;
132 double wl0, wlstep;
133 double boxpos, boxstep;
134 /* check arguments */
135 if ((nwl <= 1) | (spec == NULL) | (wlmin >= C_CMAXWL) |
136 (wlmax <= C_CMINWL) | (wlmin >= wlmax))
137 return(0.);
138 wlstep = (wlmax - wlmin)/(nwl-1.);
139 while (wlmin < C_CMINWL) {
140 wlmin += wlstep;
141 --nwl; ++spec;
142 }
143 while (wlmax > C_CMAXWL) {
144 wlmax -= wlstep;
145 --nwl;
146 }
147 if ((nwl <= 1) | (wlmin >= wlmax))
148 return(0.);
149 imax = nwl; /* box filter if necessary */
150 boxpos = 0;
151 boxstep = 1;
152 if (wlstep < C_CWLI) {
153 imax = (wlmax - wlmin)/C_CWLI + 1e-7;
154 boxpos = (wlmin - C_CMINWL)/C_CWLI;
155 boxstep = wlstep/C_CWLI;
156 wlstep = C_CWLI;
157 }
158 scale = 0.; /* get values and maximum */
159 yval = 0.;
160 pos = 0;
161 for (i = 0; i < imax; i++) {
162 va[i] = 0.; n = 0;
163 while (boxpos < i+.5 && pos < nwl) {
164 va[i] += spec[pos++];
165 n++;
166 boxpos += boxstep;
167 }
168 if (n > 1)
169 va[i] /= (double)n;
170 if (va[i] > scale)
171 scale = va[i];
172 else if (va[i] < -scale)
173 scale = -va[i];
174 yval += va[i] * c_y31.ssamp[i];
175 }
176 if (scale <= 1e-7)
177 return(0.);
178 yval /= (double)c_y31.ssum;
179 scale = C_CMAXV / scale;
180 clr->ssum = 0; /* convert to our spacing */
181 wl0 = wlmin;
182 pos = 0;
183 for (i = 0, wl = C_CMINWL; i < C_CNSS; i++, wl += C_CWLI)
184 if ((wl < wlmin) | (wl > wlmax))
185 clr->ssamp[i] = 0;
186 else {
187 while (wl0 + wlstep < wl+.01) {
188 wl0 += wlstep;
189 pos++;
190 }
191 if ((wl+.01 >= wl0) & (wl-.01 <= wl0))
192 clr->ssamp[i] = scale*va[pos] + frand();
193 else /* interpolate if necessary */
194 clr->ssamp[i] = frand() + scale / wlstep *
195 ( va[pos]*(wl0+wlstep - wl) +
196 va[pos+1]*(wl - wl0) );
197 clr->ssum += clr->ssamp[i];
198 }
199 clr->flags = C_CDSPEC|C_CSSPEC;
200 return(yval);
201 }
202
203
204 /* check if colors are equivalent */
205 int
206 c_equiv(C_COLOR *c1, C_COLOR *c2)
207 {
208 long thresh;
209 int i;
210
211 if (c1 == c2)
212 return(1);
213 if ((c1->flags ^ c2->flags) & (C_CDXY|C_CDSPEC))
214 return(0); /* defined differently, so no */
215 c_ccvt(c1, C_CSXY);
216 c_ccvt(c2, C_CSXY); /* first check chromaticities */
217 if (fabs(c1->cx - c2->cx) + fabs(c1->cy - c2->cy) > .015)
218 return(0); /* mismatch means definitely different */
219 if (c1->flags & C_CDXY)
220 return(1); /* done if defined as (x,y) */
221 thresh = C_CMAXV/200*(c1->ssum + c2->ssum);
222 for (i = 0; i < C_CNSS; i++) /* else compare spectra */
223 if (labs(c1->ssamp[i]*c2->ssum - c2->ssamp[i]*c1->ssum) > thresh)
224 return(0);
225 return(1);
226 }
227
228
229 /* check if color is grey */
230 int
231 c_isgrey(C_COLOR *clr)
232 {
233 if (!(clr->flags & (C_CSXY|C_CSSPEC)))
234 return(1); /* no settings == grey */
235
236 c_ccvt(clr, C_CSXY);
237
238 return(clr->cx >= .323 && clr->cx <= .343 &&
239 clr->cy >= .323 && clr->cy <= .343);
240 }
241
242 /* convert color representations */
243 void
244 c_ccvt(C_COLOR *clr, int fl)
245 {
246 double x, y, z;
247 int i;
248
249 fl &= ~clr->flags; /* ignore what's done */
250 if (!fl) /* everything's done! */
251 return;
252 if (!(clr->flags & (C_CSXY|C_CSSPEC))) {
253 *clr = c_dfcolor; /* nothing was set! */
254 return;
255 }
256 if (fl & C_CSXY) { /* cspec -> cxy */
257 x = y = z = 0.;
258 for (i = 0; i < C_CNSS; i++) {
259 x += c_x31.ssamp[i] * clr->ssamp[i];
260 y += c_y31.ssamp[i] * clr->ssamp[i];
261 z += c_z31.ssamp[i] * clr->ssamp[i];
262 }
263 x /= (double)c_x31.ssum;
264 y /= (double)c_y31.ssum;
265 z /= (double)c_z31.ssum;
266 z += x + y;
267 if (z > 1e-6) {
268 clr->cx = x / z;
269 clr->cy = y / z;
270 } else
271 clr->cx = clr->cy = 1./3.;
272 clr->flags |= C_CSXY;
273 }
274 if (fl & C_CSSPEC) { /* cxy -> cspec */
275 x = clr->cx;
276 y = clr->cy;
277 z = 1. - x - y;
278 clr->ssum = 0;
279 for (i = 0; i < C_CNSS; i++) {
280 clr->ssamp[i] = x*cie_xp.ssamp[i] + y*cie_yp.ssamp[i]
281 + z*cie_zp.ssamp[i] + frand();
282 if (clr->ssamp[i] < 0) /* out of gamut! */
283 clr->ssamp[i] = 0;
284 else
285 clr->ssum += clr->ssamp[i];
286 }
287 clr->flags |= C_CSSPEC;
288 }
289 if (fl & C_CSEFF) { /* compute efficacy */
290 if (clr->flags & C_CSSPEC) { /* from spectrum */
291 y = 0.;
292 for (i = 0; i < C_CNSS; i++)
293 y += c_y31.ssamp[i] * clr->ssamp[i];
294 clr->eff = C_CLPWM * y / (clr->ssum + 0.0001);
295 } else /* clr->flags & C_CSXY */ { /* from (x,y) */
296 clr->eff = clr->cx*c_x31.eff + clr->cy*c_y31.eff +
297 (1. - clr->cx - clr->cy)*c_z31.eff;
298 }
299 clr->flags |= C_CSEFF;
300 }
301 }
302
303 /* mix two colors according to weights given -- cres may be c1 or c2 */
304 void
305 c_cmix(C_COLOR *cres, double w1, C_COLOR *c1, double w2, C_COLOR *c2)
306 {
307 double scale;
308 int i;
309
310 if (w1 == 0) {
311 *cres = *c2;
312 return;
313 }
314 if (w2 == 0) {
315 *cres = *c1;
316 return;
317 }
318 if ((c1->flags|c2->flags) & C_CDSPEC) { /* spectral mixing */
319 float cmix[C_CNSS];
320
321 c_ccvt(c1, C_CSSPEC|C_CSEFF);
322 c_ccvt(c2, C_CSSPEC|C_CSEFF);
323 if (c1->ssum*c2->ssum == 0) {
324 *cres = c_dfcolor;
325 return;
326 }
327 w1 /= c1->eff*c1->ssum;
328 w2 /= c2->eff*c2->ssum;
329 scale = 0.;
330 for (i = 0; i < C_CNSS; i++) {
331 cmix[i] = w1*c1->ssamp[i] + w2*c2->ssamp[i];
332 if (cmix[i] > scale)
333 scale = cmix[i];
334 else if (cmix[i] < -scale)
335 scale = -cmix[i];
336 }
337 scale = C_CMAXV / scale;
338 cres->ssum = 0;
339 for (i = 0; i < C_CNSS; i++)
340 cres->ssum += cres->ssamp[i] = scale*cmix[i] + frand();
341 cres->flags = C_CDSPEC|C_CSSPEC;
342 } else { /* CIE xy mixing */
343 c_ccvt(c1, C_CSXY);
344 c_ccvt(c2, C_CSXY);
345 w1 *= c2->cy;
346 w2 *= c1->cy;
347 scale = w1 + w2;
348 if (scale == 0.) {
349 *cres = c_dfcolor;
350 return;
351 }
352 scale = 1. / scale;
353 cres->cx = (w1*c1->cx + w2*c2->cx) * scale;
354 cres->cy = (w1*c1->cy + w2*c2->cy) * scale;
355 cres->flags = C_CDXY|C_CSXY;
356 }
357 }
358
359 /* multiply two colors -- cres may be c1 or c2 */
360 double
361 c_cmult(C_COLOR *cres, C_COLOR *c1, double y1, C_COLOR *c2, double y2)
362 {
363 double yres;
364 int i;
365
366 if ((c1->flags|c2->flags) & C_CDSPEC) {
367 long cmix[C_CNSS], cmax = 0; /* spectral multiply */
368
369 c_ccvt(c1, C_CSSPEC|C_CSEFF);
370 c_ccvt(c2, C_CSSPEC|C_CSEFF);
371 for (i = 0; i < C_CNSS; i++) {
372 cmix[i] = c1->ssamp[i] * c2->ssamp[i];
373 if (cmix[i] > cmax)
374 cmax = cmix[i];
375 else if (cmix[i] < -cmax)
376 cmax = -cmix[i];
377 }
378 cmax /= C_CMAXV;
379 if (!cmax) {
380 *cres = c_dfcolor;
381 return(0.);
382 }
383 cres->ssum = 0;
384 for (i = 0; i < C_CNSS; i++)
385 cres->ssum += cres->ssamp[i] = cmix[i] / cmax;
386 cres->flags = C_CDSPEC|C_CSSPEC;
387
388 c_ccvt(cres, C_CSEFF); /* below is correct */
389 yres = y1 * y2 * c_y31.ssum * C_CLPWM /
390 (c1->eff*c1->ssum * c2->eff*c2->ssum) *
391 cres->eff*( cres->ssum*(double)cmax +
392 C_CNSS/2.0*(cmax-1) );
393 } else {
394 float rgb1[3], rgb2[3]; /* CIE xy multiply */
395
396 c_toSharpRGB(c1, y1, rgb1);
397 c_toSharpRGB(c2, y2, rgb2);
398
399 rgb2[0] *= rgb1[0]; rgb2[1] *= rgb1[1]; rgb2[2] *= rgb1[2];
400
401 yres = c_fromSharpRGB(rgb2, cres);
402 }
403 return(yres);
404 }
405
406 #define C1 3.741832e-16 /* W-m^2 */
407 #define C2 1.4388e-2 /* m-K */
408
409 #define bbsp(l,t) (C1/((l)*(l)*(l)*(l)*(l)*(exp(C2/((t)*(l)))-1.)))
410 #define bblm(t) (C2*0.2/(t))
411
412 /* set black body spectrum */
413 int
414 c_bbtemp(C_COLOR *clr, double tk)
415 {
416 double sf, wl;
417 int i;
418
419 if (tk < 1000)
420 return(0);
421 wl = bblm(tk); /* scalefactor based on peak */
422 if (wl < C_CMINWL*1e-9)
423 wl = C_CMINWL*1e-9;
424 else if (wl > C_CMAXWL*1e-9)
425 wl = C_CMAXWL*1e-9;
426 sf = C_CMAXV/bbsp(wl,tk);
427 clr->ssum = 0;
428 for (i = 0; i < C_CNSS; i++) {
429 wl = (C_CMINWL + i*C_CWLI)*1e-9;
430 clr->ssum += clr->ssamp[i] = sf*bbsp(wl,tk) + frand();
431 }
432 clr->flags = C_CDSPEC|C_CSSPEC;
433 return(1);
434 }
435
436 #undef C1
437 #undef C2
438 #undef bbsp
439 #undef bblm
440
441 #define UV_NORMF 410.
442
443 /* encode (x,y) chromaticity */
444 C_CHROMA
445 c_encodeChroma(C_COLOR *clr)
446 {
447 double df;
448 int ub, vb;
449
450 c_ccvt(clr, C_CSXY);
451 df = UV_NORMF/(-2.*clr->cx + 12.*clr->cy + 3.);
452 ub = 4.*clr->cx*df + frand();
453 if (ub > 0xff) ub = 0xff;
454 else ub *= (ub > 0);
455 vb = 9.*clr->cy*df + frand();
456 if (vb > 0xff) vb = 0xff;
457 else vb *= (vb > 0);
458
459 return(vb<<8 | ub);
460 }
461
462 /* decode (x,y) chromaticity */
463 void
464 c_decodeChroma(C_COLOR *cres, C_CHROMA ccode)
465 {
466 double up = (ccode & 0xff)*(1./UV_NORMF);
467 double vp = (ccode>>8 & 0xff)*(1./UV_NORMF);
468 double df = 1./(6.*up - 16.*vp + 12.);
469
470 cres->cx = 9.*up * df;
471 cres->cy = 4.*vp * df;
472 cres->flags = C_CDXY|C_CSXY;
473 }
474
475 #undef UV_NORMF