ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.11
Committed: Sat Feb 22 02:07:22 2003 UTC (21 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.10: +177 -11 lines
Log Message:
Changes and check-in for 3.5 release
Includes new source files and modifications not recorded for many years
See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id$";
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 /* ====================================================================
12 * The Radiance Software License, Version 1.0
13 *
14 * Copyright (c) 1990 - 2002 The Regents of the University of California,
15 * through Lawrence Berkeley National Laboratory. All rights reserved.
16 *
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions
19 * are met:
20 *
21 * 1. Redistributions of source code must retain the above copyright
22 * notice, this list of conditions and the following disclaimer.
23 *
24 * 2. Redistributions in binary form must reproduce the above copyright
25 * notice, this list of conditions and the following disclaimer in
26 * the documentation and/or other materials provided with the
27 * distribution.
28 *
29 * 3. The end-user documentation included with the redistribution,
30 * if any, must include the following acknowledgment:
31 * "This product includes Radiance software
32 * (http://radsite.lbl.gov/)
33 * developed by the Lawrence Berkeley National Laboratory
34 * (http://www.lbl.gov/)."
35 * Alternately, this acknowledgment may appear in the software itself,
36 * if and wherever such third-party acknowledgments normally appear.
37 *
38 * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
39 * and "The Regents of the University of California" must
40 * not be used to endorse or promote products derived from this
41 * software without prior written permission. For written
42 * permission, please contact [email protected].
43 *
44 * 5. Products derived from this software may not be called "Radiance",
45 * nor may "Radiance" appear in their name, without prior written
46 * permission of Lawrence Berkeley National Laboratory.
47 *
48 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
49 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
50 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
51 * DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
52 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
53 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
54 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
55 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
56 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
57 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
58 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
59 * SUCH DAMAGE.
60 * ====================================================================
61 *
62 * This software consists of voluntary contributions made by many
63 * individuals on behalf of Lawrence Berkeley National Laboratory. For more
64 * information on Lawrence Berkeley National Laboratory, please see
65 * <http://www.lbl.gov/>.
66 */
67
68 #include "color.h"
69 #include <string.h>
70
71 #define CEPS 1e-4 /* color epsilon */
72
73 #define CEQ(v1,v2) ((v1) <= (v2)+CEPS && (v2) <= (v1)+CEPS)
74
75 #define XYEQ(c1,c2) (CEQ((c1)[CIEX],(c2)[CIEX]) && CEQ((c1)[CIEY],(c2)[CIEY]))
76
77
78 RGBPRIMS stdprims = STDPRIMS; /* standard primary chromaticities */
79
80 COLOR cblack = BLKCOLOR; /* global black color */
81 COLOR cwhite = WHTCOLOR; /* global white color */
82
83 float xyneu[2] = {1./3., 1./3.}; /* neutral xy chromaticities */
84
85 /*
86 * The following table contains the CIE tristimulus integrals
87 * for X, Y, and Z. The table is cumulative, so that
88 * each color coordinate integrates to 1.
89 */
90
91 #define STARTWL 380 /* starting wavelength (nanometers) */
92 #define INCWL 10 /* wavelength increment */
93 #define NINC 40 /* # of values */
94
95 static BYTE chroma[3][NINC] = {
96 { /* X */
97 0, 0, 0, 2, 6, 13, 22, 30, 36, 41,
98 42, 43, 43, 44, 46, 52, 60, 71, 87, 106,
99 128, 153, 178, 200, 219, 233, 243, 249, 252, 254,
100 255, 255, 255, 255, 255, 255, 255, 255, 255, 255
101 }, { /* Y */
102 0, 0, 0, 0, 0, 1, 2, 4, 7, 11,
103 17, 24, 34, 48, 64, 84, 105, 127, 148, 169,
104 188, 205, 220, 232, 240, 246, 250, 253, 254, 255,
105 255, 255, 255, 255, 255, 255, 255, 255, 255, 255
106 }, { /* Z */
107 0, 0, 2, 10, 32, 66, 118, 153, 191, 220,
108 237, 246, 251, 253, 254, 255, 255, 255, 255, 255,
109 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
110 255, 255, 255, 255, 255, 255, 255, 255, 255, 255
111 }
112 };
113
114 COLORMAT xyz2rgbmat = { /* XYZ to RGB (no white balance) */
115 {(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g)/CIE_C_rD,
116 (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b)/CIE_C_rD,
117 (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g)/CIE_C_rD},
118 {(CIE_y_b - CIE_y_r - CIE_y_b*CIE_x_r + CIE_y_r*CIE_x_b)/CIE_C_gD,
119 (CIE_x_r - CIE_x_b - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r)/CIE_C_gD,
120 (CIE_x_b*CIE_y_r - CIE_x_r*CIE_y_b)/CIE_C_gD},
121 {(CIE_y_r - CIE_y_g - CIE_y_r*CIE_x_g + CIE_y_g*CIE_x_r)/CIE_C_bD,
122 (CIE_x_g - CIE_x_r - CIE_x_g*CIE_y_r + CIE_x_r*CIE_y_g)/CIE_C_bD,
123 (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r)/CIE_C_bD}
124 };
125
126 COLORMAT rgb2xyzmat = { /* RGB to XYZ (no white balance) */
127 {CIE_x_r*CIE_C_rD/CIE_D,CIE_x_g*CIE_C_gD/CIE_D,CIE_x_b*CIE_C_bD/CIE_D},
128 {CIE_y_r*CIE_C_rD/CIE_D,CIE_y_g*CIE_C_gD/CIE_D,CIE_y_b*CIE_C_bD/CIE_D},
129 {(1.-CIE_x_r-CIE_y_r)*CIE_C_rD/CIE_D,
130 (1.-CIE_x_g-CIE_y_g)*CIE_C_gD/CIE_D,
131 (1.-CIE_x_b-CIE_y_b)*CIE_C_bD/CIE_D}
132 };
133
134 COLORMAT vkmat = { /* Sharp primary matrix */
135 { 1.2694, -0.0988, -0.1706},
136 {-0.8364, 1.8006, 0.0357},
137 { 0.0297, -0.0315, 1.0018}
138 };
139
140 COLORMAT ivkmat = { /* inverse Sharp primary matrix */
141 { 0.8156, 0.0472, 0.1372},
142 { 0.3791, 0.5769, 0.0440},
143 {-0.0123, 0.0167, 0.9955}
144 };
145
146
147 void
148 spec_rgb(col, s, e) /* compute RGB color from spectral range */
149 COLOR col;
150 int s, e;
151 {
152 COLOR ciecolor;
153
154 spec_cie(ciecolor, s, e);
155 cie_rgb(col, ciecolor);
156 }
157
158
159 void
160 spec_cie(col, s, e) /* compute a color from a spectral range */
161 COLOR col; /* returned color */
162 int s, e; /* starting and ending wavelengths */
163 {
164 register int i, d, r;
165
166 s -= STARTWL;
167 if (s < 0)
168 s = 0;
169
170 e -= STARTWL;
171 if (e <= s) {
172 col[CIEX] = col[CIEY] = col[CIEZ] = 0.0;
173 return;
174 }
175 if (e >= INCWL*(NINC - 1))
176 e = INCWL*(NINC - 1) - 1;
177
178 d = e / INCWL; /* interpolate values */
179 r = e % INCWL;
180 for (i = 0; i < 3; i++)
181 col[i] = chroma[i][d]*(INCWL - r) + chroma[i][d + 1]*r;
182
183 d = s / INCWL;
184 r = s % INCWL;
185 for (i = 0; i < 3; i++)
186 col[i] -= chroma[i][d]*(INCWL - r) - chroma[i][d + 1]*r;
187
188 col[CIEX] = (col[CIEX] + 0.5) * (1./(256*INCWL));
189 col[CIEY] = (col[CIEY] + 0.5) * (1./(256*INCWL));
190 col[CIEZ] = (col[CIEZ] + 0.5) * (1./(256*INCWL));
191 }
192
193
194 void
195 cie_rgb(rgb, xyz) /* convert CIE color to standard RGB */
196 COLOR rgb;
197 COLOR xyz;
198 {
199 colortrans(rgb, xyz2rgbmat, xyz);
200 clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
201 }
202
203
204 int
205 clipgamut(col, brt, gamut, lower, upper) /* clip to gamut cube */
206 COLOR col;
207 double brt;
208 int gamut;
209 COLOR lower, upper;
210 {
211 int rflags = 0;
212 double brtmin, brtmax, v, vv;
213 COLOR cgry;
214 register int i;
215 /* check for no check */
216 if (gamut == 0) return(0);
217 /* check brightness limits */
218 brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
219 if (gamut & CGAMUT_LOWER && brt < brtmin) {
220 copycolor(col, lower);
221 return(CGAMUT_LOWER);
222 }
223 brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
224 if (gamut & CGAMUT_UPPER && brt > brtmax) {
225 copycolor(col, upper);
226 return(CGAMUT_UPPER);
227 }
228 /* compute equivalent grey */
229 v = (brt - brtmin)/(brtmax - brtmin);
230 for (i = 0; i < 3; i++)
231 cgry[i] = v*upper[i] + (1.-v)*lower[i];
232 vv = 1.; /* check each limit */
233 for (i = 0; i < 3; i++)
234 if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
235 v = (lower[i]+CEPS - cgry[i])/(col[i] - cgry[i]);
236 if (v < vv) vv = v;
237 rflags |= CGAMUT_LOWER;
238 } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
239 v = (upper[i]-CEPS - cgry[i])/(col[i] - cgry[i]);
240 if (v < vv) vv = v;
241 rflags |= CGAMUT_UPPER;
242 }
243 if (rflags) /* desaturate to cube face */
244 for (i = 0; i < 3; i++)
245 col[i] = vv*col[i] + (1.-vv)*cgry[i];
246 return(rflags);
247 }
248
249
250 void
251 colortrans(c2, mat, c1) /* convert c1 by mat and put into c2 */
252 register COLOR c2;
253 register COLORMAT mat;
254 register COLOR c1;
255 {
256 COLOR cout;
257
258 cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
259 cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
260 cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];
261
262 copycolor(c2, cout);
263 }
264
265
266 void
267 multcolormat(m3, m2, m1) /* multiply m1 by m2 and put into m3 */
268 COLORMAT m3; /* m3 can be either m1 or m2 w/o harm */
269 COLORMAT m2, m1;
270 {
271 COLORMAT mt;
272 register int i, j;
273
274 for (i = 0; i < 3; i++)
275 for (j = 0; j < 3; j++)
276 mt[i][j] = m1[i][0]*m2[0][j] +
277 m1[i][1]*m2[1][j] +
278 m1[i][2]*m2[2][j] ;
279 cpcolormat(m3, mt);
280 }
281
282
283 void
284 compxyz2rgbmat(mat, pr) /* compute conversion from CIE to RGB space */
285 COLORMAT mat;
286 register RGBPRIMS pr;
287 {
288 double C_rD, C_gD, C_bD;
289
290 if (pr == stdprims) { /* can use xyz2rgbmat */
291 cpcolormat(mat, xyz2rgbmat);
292 return;
293 }
294 C_rD = (1./pr[WHT][CIEY]) *
295 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
296 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
297 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
298 C_gD = (1./pr[WHT][CIEY]) *
299 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
300 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
301 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
302 C_bD = (1./pr[WHT][CIEY]) *
303 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
304 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
305 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
306
307 mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
308 pr[BLU][CIEX]*pr[GRN][CIEY] +
309 pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
310 mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
311 pr[BLU][CIEX]*pr[GRN][CIEY] +
312 pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
313 mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
314 pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
315 mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
316 pr[BLU][CIEY]*pr[RED][CIEX] +
317 pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
318 mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
319 pr[RED][CIEX]*pr[BLU][CIEY] +
320 pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
321 mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
322 pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
323 mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
324 pr[RED][CIEY]*pr[GRN][CIEX] +
325 pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
326 mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
327 pr[GRN][CIEX]*pr[RED][CIEY] +
328 pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
329 mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
330 pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
331 }
332
333
334 void
335 comprgb2xyzmat(mat, pr) /* compute conversion from RGB to CIE space */
336 COLORMAT mat;
337 register RGBPRIMS pr;
338 {
339 double C_rD, C_gD, C_bD, D;
340
341 if (pr == stdprims) { /* can use rgb2xyzmat */
342 cpcolormat(mat, rgb2xyzmat);
343 return;
344 }
345 C_rD = (1./pr[WHT][CIEY]) *
346 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
347 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
348 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
349 C_gD = (1./pr[WHT][CIEY]) *
350 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
351 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
352 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
353 C_bD = (1./pr[WHT][CIEY]) *
354 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
355 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
356 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
357 D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
358 pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
359 pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
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 }
370
371
372 void
373 comprgb2rgbmat(mat, pr1, pr2) /* compute conversion from RGB1 to RGB2 */
374 COLORMAT mat;
375 RGBPRIMS pr1, pr2;
376 {
377 COLORMAT pr1toxyz, xyztopr2;
378
379 if (pr1 == pr2) {
380 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
381 mat[0][1] = mat[0][2] = mat[1][0] =
382 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
383 return;
384 }
385 comprgb2xyzmat(pr1toxyz, pr1);
386 compxyz2rgbmat(xyztopr2, pr2);
387 /* combine transforms */
388 multcolormat(mat, pr1toxyz, xyztopr2);
389 }
390
391
392 void
393 compxyzWBmat(mat, wht1, wht2) /* CIE von Kries transform from wht1 to wht2 */
394 COLORMAT mat;
395 float wht1[2], wht2[2];
396 {
397 COLOR cw1, cw2;
398 if (XYEQ(wht1,wht2)) {
399 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
400 mat[0][1] = mat[0][2] = mat[1][0] =
401 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
402 return;
403 }
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 cw2[RED] = wht2[CIEX]/wht2[CIEY];
409 cw2[GRN] = 1.;
410 cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
411 colortrans(cw2, vkmat, cw2);
412 mat[0][0] = cw2[RED]/cw1[RED];
413 mat[1][1] = cw2[GRN]/cw1[GRN];
414 mat[2][2] = cw2[BLU]/cw1[BLU];
415 mat[0][1] = mat[0][2] = mat[1][0] =
416 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
417 multcolormat(mat, vkmat, mat);
418 multcolormat(mat, mat, ivkmat);
419 }
420
421
422 void
423 compxyz2rgbWBmat(mat, pr) /* von Kries conversion from CIE to RGB space */
424 COLORMAT mat;
425 RGBPRIMS pr;
426 {
427 COLORMAT wbmat;
428
429 compxyz2rgbmat(mat, pr);
430 if (XYEQ(pr[WHT],xyneu))
431 return;
432 compxyzWBmat(wbmat, xyneu, pr[WHT]);
433 multcolormat(mat, wbmat, mat);
434 }
435
436 void
437 comprgb2xyzWBmat(mat, pr) /* von Kries conversion from RGB to CIE space */
438 COLORMAT mat;
439 RGBPRIMS pr;
440 {
441 COLORMAT wbmat;
442
443 comprgb2xyzmat(mat, pr);
444 if (XYEQ(pr[WHT],xyneu))
445 return;
446 compxyzWBmat(wbmat, pr[WHT], xyneu);
447 multcolormat(mat, mat, wbmat);
448 }
449
450 void
451 comprgb2rgbWBmat(mat, pr1, pr2) /* von Kries conversion from RGB1 to RGB2 */
452 COLORMAT mat;
453 RGBPRIMS pr1, pr2;
454 {
455 COLORMAT pr1toxyz, xyztopr2, wbmat;
456
457 if (pr1 == pr2) {
458 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
459 mat[0][1] = mat[0][2] = mat[1][0] =
460 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
461 return;
462 }
463 comprgb2xyzmat(pr1toxyz, pr1);
464 compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]);
465 compxyz2rgbmat(xyztopr2, pr2);
466 /* combine transforms */
467 multcolormat(mat, pr1toxyz, wbmat);
468 multcolormat(mat, mat, xyztopr2);
469 }