ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.17
Committed: Mon Nov 2 23:45:11 2009 UTC (14 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.16: +5 -3 lines
Log Message:
Made color primitives check a little more forgiving

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: spec_rgb.c,v 2.16 2007/10/08 18:07:55 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;
245
246 for (i = 0; i < 4; i++) {
247 if ((pr[i][CIEX] <= -CEPS) | (pr[i][CIEY] <= -CEPS))
248 return(0);
249 if ((pr[i][CIEX] >= 1.+CEPS) | (pr[i][CIEY] >= 1.+CEPS))
250 return(0);
251 if (pr[i][CIEX] + pr[i][CIEY] >= 1.+CEPS)
252 return(0);
253 }
254 return(1);
255 }
256
257
258
259 int
260 compxyz2rgbmat( /* compute conversion from CIE to RGB space */
261 COLORMAT mat,
262 register RGBPRIMS pr
263 )
264 {
265 double C_rD, C_gD, C_bD;
266
267 if (pr == stdprims) { /* can use xyz2rgbmat */
268 cpcolormat(mat, xyz2rgbmat);
269 return(1);
270 }
271 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
272 return(0);
273 C_rD = (1./pr[WHT][CIEY]) *
274 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
275 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
276 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
277 if (CEQ(C_rD,0.))
278 return(0);
279 C_gD = (1./pr[WHT][CIEY]) *
280 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
281 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
282 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
283 if (CEQ(C_gD,0.))
284 return(0);
285 C_bD = (1./pr[WHT][CIEY]) *
286 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
287 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
288 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
289 if (CEQ(C_bD,0.))
290 return(0);
291 mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
292 pr[BLU][CIEX]*pr[GRN][CIEY] +
293 pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
294 mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
295 pr[BLU][CIEX]*pr[GRN][CIEY] +
296 pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
297 mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
298 pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
299 mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
300 pr[BLU][CIEY]*pr[RED][CIEX] +
301 pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
302 mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
303 pr[RED][CIEX]*pr[BLU][CIEY] +
304 pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
305 mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
306 pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
307 mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
308 pr[RED][CIEY]*pr[GRN][CIEX] +
309 pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
310 mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
311 pr[GRN][CIEX]*pr[RED][CIEY] +
312 pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
313 mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
314 pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
315 return(1);
316 }
317
318
319 int
320 comprgb2xyzmat( /* compute conversion from RGB to CIE space */
321 COLORMAT mat,
322 register RGBPRIMS pr
323 )
324 {
325 double C_rD, C_gD, C_bD, D;
326
327 if (pr == stdprims) { /* can use rgb2xyzmat */
328 cpcolormat(mat, rgb2xyzmat);
329 return(1);
330 }
331 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
332 return(0);
333 C_rD = (1./pr[WHT][CIEY]) *
334 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
335 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
336 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
337 C_gD = (1./pr[WHT][CIEY]) *
338 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
339 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
340 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
341 C_bD = (1./pr[WHT][CIEY]) *
342 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
343 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
344 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
345 D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
346 pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
347 pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
348 if (CEQ(D,0.))
349 return(0);
350 mat[0][0] = pr[RED][CIEX]*C_rD/D;
351 mat[0][1] = pr[GRN][CIEX]*C_gD/D;
352 mat[0][2] = pr[BLU][CIEX]*C_bD/D;
353 mat[1][0] = pr[RED][CIEY]*C_rD/D;
354 mat[1][1] = pr[GRN][CIEY]*C_gD/D;
355 mat[1][2] = pr[BLU][CIEY]*C_bD/D;
356 mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
357 mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
358 mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
359 return(1);
360 }
361
362
363 int
364 comprgb2rgbmat( /* compute conversion from RGB1 to RGB2 */
365 COLORMAT mat,
366 RGBPRIMS pr1,
367 RGBPRIMS pr2
368 )
369 {
370 COLORMAT pr1toxyz, xyztopr2;
371
372 if (pr1 == pr2) {
373 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
374 mat[0][1] = mat[0][2] = mat[1][0] =
375 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
376 return(1);
377 }
378 if (!comprgb2xyzmat(pr1toxyz, pr1))
379 return(0);
380 if (!compxyz2rgbmat(xyztopr2, pr2))
381 return(0);
382 /* combine transforms */
383 multcolormat(mat, pr1toxyz, xyztopr2);
384 return(1);
385 }
386
387
388 int
389 compxyzWBmat( /* CIE von Kries transform from wht1 to wht2 */
390 COLORMAT mat,
391 float wht1[2],
392 float wht2[2]
393 )
394 {
395 COLOR cw1, cw2;
396 if (XYEQ(wht1,wht2)) {
397 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
398 mat[0][1] = mat[0][2] = mat[1][0] =
399 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
400 return(1);
401 }
402 if (CEQ(wht1[CIEX],0.) | CEQ(wht1[CIEY],0.))
403 return(0);
404 cw1[RED] = wht1[CIEX]/wht1[CIEY];
405 cw1[GRN] = 1.;
406 cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
407 colortrans(cw1, vkmat, cw1);
408 if (CEQ(wht2[CIEX],0.) | CEQ(wht2[CIEY],0.))
409 return(0);
410 cw2[RED] = wht2[CIEX]/wht2[CIEY];
411 cw2[GRN] = 1.;
412 cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
413 colortrans(cw2, vkmat, cw2);
414 if (CEQ(cw1[RED],0.) | CEQ(cw1[GRN],0.) | CEQ(cw1[BLU],0.))
415 return(0);
416 mat[0][0] = cw2[RED]/cw1[RED];
417 mat[1][1] = cw2[GRN]/cw1[GRN];
418 mat[2][2] = cw2[BLU]/cw1[BLU];
419 mat[0][1] = mat[0][2] = mat[1][0] =
420 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
421 multcolormat(mat, vkmat, mat);
422 multcolormat(mat, mat, ivkmat);
423 return(1);
424 }
425
426
427 int
428 compxyz2rgbWBmat( /* von Kries conversion from CIE to RGB space */
429 COLORMAT mat,
430 RGBPRIMS pr
431 )
432 {
433 COLORMAT wbmat;
434
435 if (!compxyz2rgbmat(mat, pr))
436 return(0);
437 if (XYEQ(pr[WHT],xyneu))
438 return(1);
439 if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
440 return(0);
441 multcolormat(mat, wbmat, mat);
442 return(1);
443 }
444
445 int
446 comprgb2xyzWBmat( /* von Kries conversion from RGB to CIE space */
447 COLORMAT mat,
448 RGBPRIMS pr
449 )
450 {
451 COLORMAT wbmat;
452
453 if (!comprgb2xyzmat(mat, pr))
454 return(0);
455 if (XYEQ(pr[WHT],xyneu))
456 return(1);
457 if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
458 return(0);
459 multcolormat(mat, mat, wbmat);
460 return(1);
461 }
462
463 int
464 comprgb2rgbWBmat( /* von Kries conversion from RGB1 to RGB2 */
465 COLORMAT mat,
466 RGBPRIMS pr1,
467 RGBPRIMS pr2
468 )
469 {
470 COLORMAT pr1toxyz, xyztopr2, wbmat;
471
472 if (pr1 == pr2) {
473 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
474 mat[0][1] = mat[0][2] = mat[1][0] =
475 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
476 return(1);
477 }
478 if (!comprgb2xyzmat(pr1toxyz, pr1))
479 return(0);
480 if (!compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]))
481 return(0);
482 if (!compxyz2rgbmat(xyztopr2, pr2))
483 return(0);
484 /* combine transforms */
485 multcolormat(mat, pr1toxyz, wbmat);
486 multcolormat(mat, mat, xyztopr2);
487 return(1);
488 }