ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.24
Committed: Fri May 18 22:48:39 2012 UTC (11 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2, rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R2P1
Changes since 2.23: +5 -5 lines
Log Message:
Even wider acceptance of whacky color primaries

File Contents

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