ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.12
Committed: Tue Feb 25 02:47:22 2003 UTC (21 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R5
Changes since 2.11: +1 -56 lines
Log Message:
Replaced inline copyright notice with #include "copyright.h"

File Contents

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.11 static const char RCSid[] = "$Id$";
3 greg 1.1 #endif
4     /*
5     * Convert colors and spectral ranges.
6 greg 2.11 * Added von Kries white-balance calculations 10/01 (GW).
7     *
8     * Externals declared in color.h
9     */
10    
11 greg 2.12 #include "copyright.h"
12 greg 1.1
13     #include "color.h"
14 greg 2.11 #include <string.h>
15    
16     #define CEPS 1e-4 /* color epsilon */
17    
18     #define CEQ(v1,v2) ((v1) <= (v2)+CEPS && (v2) <= (v1)+CEPS)
19 greg 1.1
20 greg 2.11 #define XYEQ(c1,c2) (CEQ((c1)[CIEX],(c2)[CIEX]) && CEQ((c1)[CIEY],(c2)[CIEY]))
21 greg 2.5
22 gwlarson 2.10
23 greg 2.11 RGBPRIMS stdprims = STDPRIMS; /* standard primary chromaticities */
24 greg 2.5
25 greg 2.8 COLOR cblack = BLKCOLOR; /* global black color */
26     COLOR cwhite = WHTCOLOR; /* global white color */
27    
28 greg 2.11 float xyneu[2] = {1./3., 1./3.}; /* neutral xy chromaticities */
29    
30 greg 1.1 /*
31     * The following table contains the CIE tristimulus integrals
32     * for X, Y, and Z. The table is cumulative, so that
33     * each color coordinate integrates to 1.
34     */
35    
36     #define STARTWL 380 /* starting wavelength (nanometers) */
37     #define INCWL 10 /* wavelength increment */
38     #define NINC 40 /* # of values */
39    
40     static BYTE chroma[3][NINC] = {
41     { /* X */
42     0, 0, 0, 2, 6, 13, 22, 30, 36, 41,
43     42, 43, 43, 44, 46, 52, 60, 71, 87, 106,
44     128, 153, 178, 200, 219, 233, 243, 249, 252, 254,
45     255, 255, 255, 255, 255, 255, 255, 255, 255, 255
46     }, { /* Y */
47     0, 0, 0, 0, 0, 1, 2, 4, 7, 11,
48     17, 24, 34, 48, 64, 84, 105, 127, 148, 169,
49     188, 205, 220, 232, 240, 246, 250, 253, 254, 255,
50     255, 255, 255, 255, 255, 255, 255, 255, 255, 255
51     }, { /* Z */
52     0, 0, 2, 10, 32, 66, 118, 153, 191, 220,
53     237, 246, 251, 253, 254, 255, 255, 255, 255, 255,
54     255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
55     255, 255, 255, 255, 255, 255, 255, 255, 255, 255
56     }
57     };
58    
59 greg 2.11 COLORMAT xyz2rgbmat = { /* XYZ to RGB (no white balance) */
60 greg 2.3 {(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g)/CIE_C_rD,
61     (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b)/CIE_C_rD,
62     (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g)/CIE_C_rD},
63     {(CIE_y_b - CIE_y_r - CIE_y_b*CIE_x_r + CIE_y_r*CIE_x_b)/CIE_C_gD,
64     (CIE_x_r - CIE_x_b - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r)/CIE_C_gD,
65     (CIE_x_b*CIE_y_r - CIE_x_r*CIE_y_b)/CIE_C_gD},
66     {(CIE_y_r - CIE_y_g - CIE_y_r*CIE_x_g + CIE_y_g*CIE_x_r)/CIE_C_bD,
67     (CIE_x_g - CIE_x_r - CIE_x_g*CIE_y_r + CIE_x_r*CIE_y_g)/CIE_C_bD,
68     (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r)/CIE_C_bD}
69 greg 1.2 };
70 greg 1.1
71 greg 2.11 COLORMAT rgb2xyzmat = { /* RGB to XYZ (no white balance) */
72 greg 2.4 {CIE_x_r*CIE_C_rD/CIE_D,CIE_x_g*CIE_C_gD/CIE_D,CIE_x_b*CIE_C_bD/CIE_D},
73     {CIE_y_r*CIE_C_rD/CIE_D,CIE_y_g*CIE_C_gD/CIE_D,CIE_y_b*CIE_C_bD/CIE_D},
74     {(1.-CIE_x_r-CIE_y_r)*CIE_C_rD/CIE_D,
75     (1.-CIE_x_g-CIE_y_g)*CIE_C_gD/CIE_D,
76     (1.-CIE_x_b-CIE_y_b)*CIE_C_bD/CIE_D}
77     };
78 greg 1.2
79 greg 2.11 COLORMAT vkmat = { /* Sharp primary matrix */
80     { 1.2694, -0.0988, -0.1706},
81     {-0.8364, 1.8006, 0.0357},
82     { 0.0297, -0.0315, 1.0018}
83     };
84    
85     COLORMAT ivkmat = { /* inverse Sharp primary matrix */
86     { 0.8156, 0.0472, 0.1372},
87     { 0.3791, 0.5769, 0.0440},
88     {-0.0123, 0.0167, 0.9955}
89     };
90 greg 1.2
91 greg 2.4
92 greg 2.11 void
93 greg 1.1 spec_rgb(col, s, e) /* compute RGB color from spectral range */
94     COLOR col;
95     int s, e;
96     {
97     COLOR ciecolor;
98    
99     spec_cie(ciecolor, s, e);
100     cie_rgb(col, ciecolor);
101     }
102    
103    
104 greg 2.11 void
105 greg 1.1 spec_cie(col, s, e) /* compute a color from a spectral range */
106     COLOR col; /* returned color */
107     int s, e; /* starting and ending wavelengths */
108     {
109     register int i, d, r;
110    
111     s -= STARTWL;
112     if (s < 0)
113     s = 0;
114    
115     e -= STARTWL;
116 greg 1.2 if (e <= s) {
117 greg 2.5 col[CIEX] = col[CIEY] = col[CIEZ] = 0.0;
118 greg 1.2 return;
119     }
120 greg 1.1 if (e >= INCWL*(NINC - 1))
121     e = INCWL*(NINC - 1) - 1;
122    
123     d = e / INCWL; /* interpolate values */
124     r = e % INCWL;
125     for (i = 0; i < 3; i++)
126     col[i] = chroma[i][d]*(INCWL - r) + chroma[i][d + 1]*r;
127    
128     d = s / INCWL;
129     r = s % INCWL;
130     for (i = 0; i < 3; i++)
131     col[i] -= chroma[i][d]*(INCWL - r) - chroma[i][d + 1]*r;
132    
133 greg 2.5 col[CIEX] = (col[CIEX] + 0.5) * (1./(256*INCWL));
134     col[CIEY] = (col[CIEY] + 0.5) * (1./(256*INCWL));
135     col[CIEZ] = (col[CIEZ] + 0.5) * (1./(256*INCWL));
136 greg 1.1 }
137    
138    
139 greg 2.11 void
140 greg 2.8 cie_rgb(rgb, xyz) /* convert CIE color to standard RGB */
141 greg 2.11 COLOR rgb;
142     COLOR xyz;
143 greg 2.8 {
144     colortrans(rgb, xyz2rgbmat, xyz);
145     clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
146     }
147    
148    
149     int
150     clipgamut(col, brt, gamut, lower, upper) /* clip to gamut cube */
151     COLOR col;
152     double brt;
153     int gamut;
154     COLOR lower, upper;
155     {
156     int rflags = 0;
157     double brtmin, brtmax, v, vv;
158     COLOR cgry;
159     register int i;
160     /* check for no check */
161     if (gamut == 0) return(0);
162     /* check brightness limits */
163     brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
164     if (gamut & CGAMUT_LOWER && brt < brtmin) {
165     copycolor(col, lower);
166     return(CGAMUT_LOWER);
167     }
168     brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
169     if (gamut & CGAMUT_UPPER && brt > brtmax) {
170     copycolor(col, upper);
171     return(CGAMUT_UPPER);
172     }
173     /* compute equivalent grey */
174     v = (brt - brtmin)/(brtmax - brtmin);
175     for (i = 0; i < 3; i++)
176     cgry[i] = v*upper[i] + (1.-v)*lower[i];
177     vv = 1.; /* check each limit */
178     for (i = 0; i < 3; i++)
179     if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
180 gwlarson 2.10 v = (lower[i]+CEPS - cgry[i])/(col[i] - cgry[i]);
181 greg 2.8 if (v < vv) vv = v;
182     rflags |= CGAMUT_LOWER;
183     } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
184 gwlarson 2.10 v = (upper[i]-CEPS - cgry[i])/(col[i] - cgry[i]);
185 greg 2.8 if (v < vv) vv = v;
186     rflags |= CGAMUT_UPPER;
187     }
188     if (rflags) /* desaturate to cube face */
189     for (i = 0; i < 3; i++)
190     col[i] = vv*col[i] + (1.-vv)*cgry[i];
191     return(rflags);
192     }
193    
194    
195 greg 2.11 void
196 greg 2.8 colortrans(c2, mat, c1) /* convert c1 by mat and put into c2 */
197 greg 2.11 register COLOR c2;
198 greg 2.5 register COLORMAT mat;
199 greg 2.11 register COLOR c1;
200 greg 1.1 {
201 greg 2.9 COLOR cout;
202    
203     cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
204     cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
205     cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];
206    
207     copycolor(c2, cout);
208 greg 2.4 }
209    
210    
211 greg 2.11 void
212 greg 2.5 multcolormat(m3, m2, m1) /* multiply m1 by m2 and put into m3 */
213 greg 2.11 COLORMAT m3; /* m3 can be either m1 or m2 w/o harm */
214     COLORMAT m2, m1;
215 greg 2.4 {
216 greg 2.5 COLORMAT mt;
217     register int i, j;
218 greg 2.4
219     for (i = 0; i < 3; i++)
220 greg 2.5 for (j = 0; j < 3; j++)
221     mt[i][j] = m1[i][0]*m2[0][j] +
222     m1[i][1]*m2[1][j] +
223     m1[i][2]*m2[2][j] ;
224     cpcolormat(m3, mt);
225     }
226    
227    
228 greg 2.11 void
229 greg 2.5 compxyz2rgbmat(mat, pr) /* compute conversion from CIE to RGB space */
230     COLORMAT mat;
231     register RGBPRIMS pr;
232     {
233     double C_rD, C_gD, C_bD;
234    
235 greg 2.6 if (pr == stdprims) { /* can use xyz2rgbmat */
236     cpcolormat(mat, xyz2rgbmat);
237     return;
238     }
239 greg 2.5 C_rD = (1./pr[WHT][CIEY]) *
240     ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
241     pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
242     pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
243     C_gD = (1./pr[WHT][CIEY]) *
244     ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
245     pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
246     pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
247     C_bD = (1./pr[WHT][CIEY]) *
248     ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
249     pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
250     pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
251    
252     mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
253     pr[BLU][CIEX]*pr[GRN][CIEY] +
254     pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
255     mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
256     pr[BLU][CIEX]*pr[GRN][CIEY] +
257     pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
258     mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
259     pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
260     mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
261     pr[BLU][CIEY]*pr[RED][CIEX] +
262     pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
263     mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
264     pr[RED][CIEX]*pr[BLU][CIEY] +
265     pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
266     mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
267     pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
268     mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
269     pr[RED][CIEY]*pr[GRN][CIEX] +
270     pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
271     mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
272     pr[GRN][CIEX]*pr[RED][CIEY] +
273     pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
274     mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
275     pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
276     }
277    
278    
279 greg 2.11 void
280 greg 2.5 comprgb2xyzmat(mat, pr) /* compute conversion from RGB to CIE space */
281     COLORMAT mat;
282     register RGBPRIMS pr;
283     {
284     double C_rD, C_gD, C_bD, D;
285    
286 greg 2.6 if (pr == stdprims) { /* can use rgb2xyzmat */
287     cpcolormat(mat, rgb2xyzmat);
288     return;
289     }
290 greg 2.5 C_rD = (1./pr[WHT][CIEY]) *
291     ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
292     pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
293     pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
294     C_gD = (1./pr[WHT][CIEY]) *
295     ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
296     pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
297     pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
298     C_bD = (1./pr[WHT][CIEY]) *
299     ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
300     pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
301     pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
302     D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
303     pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
304     pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
305     mat[0][0] = pr[RED][CIEX]*C_rD/D;
306     mat[0][1] = pr[GRN][CIEX]*C_gD/D;
307     mat[0][2] = pr[BLU][CIEX]*C_bD/D;
308     mat[1][0] = pr[RED][CIEY]*C_rD/D;
309     mat[1][1] = pr[GRN][CIEY]*C_gD/D;
310     mat[1][2] = pr[BLU][CIEY]*C_bD/D;
311     mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
312     mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
313     mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
314     }
315    
316    
317 greg 2.11 void
318 greg 2.5 comprgb2rgbmat(mat, pr1, pr2) /* compute conversion from RGB1 to RGB2 */
319     COLORMAT mat;
320     RGBPRIMS pr1, pr2;
321     {
322     COLORMAT pr1toxyz, xyztopr2;
323    
324 greg 2.6 if (pr1 == pr2) {
325     mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
326     mat[0][1] = mat[0][2] = mat[1][0] =
327     mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
328     return;
329     }
330     comprgb2xyzmat(pr1toxyz, pr1);
331     compxyz2rgbmat(xyztopr2, pr2);
332 greg 2.5 /* combine transforms */
333     multcolormat(mat, pr1toxyz, xyztopr2);
334 greg 2.11 }
335    
336    
337     void
338     compxyzWBmat(mat, wht1, wht2) /* CIE von Kries transform from wht1 to wht2 */
339     COLORMAT mat;
340     float wht1[2], wht2[2];
341     {
342     COLOR cw1, cw2;
343     if (XYEQ(wht1,wht2)) {
344     mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
345     mat[0][1] = mat[0][2] = mat[1][0] =
346     mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
347     return;
348     }
349     cw1[RED] = wht1[CIEX]/wht1[CIEY];
350     cw1[GRN] = 1.;
351     cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
352     colortrans(cw1, vkmat, cw1);
353     cw2[RED] = wht2[CIEX]/wht2[CIEY];
354     cw2[GRN] = 1.;
355     cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
356     colortrans(cw2, vkmat, cw2);
357     mat[0][0] = cw2[RED]/cw1[RED];
358     mat[1][1] = cw2[GRN]/cw1[GRN];
359     mat[2][2] = cw2[BLU]/cw1[BLU];
360     mat[0][1] = mat[0][2] = mat[1][0] =
361     mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
362     multcolormat(mat, vkmat, mat);
363     multcolormat(mat, mat, ivkmat);
364     }
365    
366    
367     void
368     compxyz2rgbWBmat(mat, pr) /* von Kries conversion from CIE to RGB space */
369     COLORMAT mat;
370     RGBPRIMS pr;
371     {
372     COLORMAT wbmat;
373    
374     compxyz2rgbmat(mat, pr);
375     if (XYEQ(pr[WHT],xyneu))
376     return;
377     compxyzWBmat(wbmat, xyneu, pr[WHT]);
378     multcolormat(mat, wbmat, mat);
379     }
380    
381     void
382     comprgb2xyzWBmat(mat, pr) /* von Kries conversion from RGB to CIE space */
383     COLORMAT mat;
384     RGBPRIMS pr;
385     {
386     COLORMAT wbmat;
387    
388     comprgb2xyzmat(mat, pr);
389     if (XYEQ(pr[WHT],xyneu))
390     return;
391     compxyzWBmat(wbmat, pr[WHT], xyneu);
392     multcolormat(mat, mat, wbmat);
393     }
394    
395     void
396     comprgb2rgbWBmat(mat, pr1, pr2) /* von Kries conversion from RGB1 to RGB2 */
397     COLORMAT mat;
398     RGBPRIMS pr1, pr2;
399     {
400     COLORMAT pr1toxyz, xyztopr2, wbmat;
401    
402     if (pr1 == pr2) {
403     mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
404     mat[0][1] = mat[0][2] = mat[1][0] =
405     mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
406     return;
407     }
408     comprgb2xyzmat(pr1toxyz, pr1);
409     compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]);
410     compxyz2rgbmat(xyztopr2, pr2);
411     /* combine transforms */
412     multcolormat(mat, pr1toxyz, wbmat);
413     multcolormat(mat, mat, xyztopr2);
414 greg 1.1 }