ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.20
Committed: Mon Jan 4 19:40:14 2010 UTC (14 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R0
Changes since 2.19: +7 -2 lines
Log Message:
Improved chromaticity check

File Contents

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