ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/ccolor.c
Revision: 3.3
Committed: Thu May 17 17:10:23 2012 UTC (11 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.2: +82 -2 lines
Log Message:
Added spectrum assignment function c_sset()

File Contents

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