ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.13
Committed: Mon Jun 30 19:04:29 2003 UTC (20 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad3R6, rad3R6P1
Changes since 2.12: +3 -2 lines
Log Message:
Removed stdio.h includes from calcomp.h, resolu.h, and color.h

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: spec_rgb.c,v 2.12 2003/02/25 02:47:22 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(col, s, e) /* compute RGB color from spectral range */
95 COLOR col;
96 int s, e;
97 {
98 COLOR ciecolor;
99
100 spec_cie(ciecolor, s, e);
101 cie_rgb(col, ciecolor);
102 }
103
104
105 void
106 spec_cie(col, s, e) /* compute a color from a spectral range */
107 COLOR col; /* returned color */
108 int s, e; /* starting and ending wavelengths */
109 {
110 register int i, d, r;
111
112 s -= STARTWL;
113 if (s < 0)
114 s = 0;
115
116 e -= STARTWL;
117 if (e <= s) {
118 col[CIEX] = col[CIEY] = col[CIEZ] = 0.0;
119 return;
120 }
121 if (e >= INCWL*(NINC - 1))
122 e = INCWL*(NINC - 1) - 1;
123
124 d = e / INCWL; /* interpolate values */
125 r = e % INCWL;
126 for (i = 0; i < 3; i++)
127 col[i] = chroma[i][d]*(INCWL - r) + chroma[i][d + 1]*r;
128
129 d = s / INCWL;
130 r = s % INCWL;
131 for (i = 0; i < 3; i++)
132 col[i] -= chroma[i][d]*(INCWL - r) - chroma[i][d + 1]*r;
133
134 col[CIEX] = (col[CIEX] + 0.5) * (1./(256*INCWL));
135 col[CIEY] = (col[CIEY] + 0.5) * (1./(256*INCWL));
136 col[CIEZ] = (col[CIEZ] + 0.5) * (1./(256*INCWL));
137 }
138
139
140 void
141 cie_rgb(rgb, xyz) /* convert CIE color to standard RGB */
142 COLOR rgb;
143 COLOR xyz;
144 {
145 colortrans(rgb, xyz2rgbmat, xyz);
146 clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
147 }
148
149
150 int
151 clipgamut(col, brt, gamut, lower, upper) /* clip to gamut cube */
152 COLOR col;
153 double brt;
154 int gamut;
155 COLOR lower, upper;
156 {
157 int rflags = 0;
158 double brtmin, brtmax, v, vv;
159 COLOR cgry;
160 register int i;
161 /* check for no check */
162 if (gamut == 0) return(0);
163 /* check brightness limits */
164 brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
165 if (gamut & CGAMUT_LOWER && brt < brtmin) {
166 copycolor(col, lower);
167 return(CGAMUT_LOWER);
168 }
169 brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
170 if (gamut & CGAMUT_UPPER && brt > brtmax) {
171 copycolor(col, upper);
172 return(CGAMUT_UPPER);
173 }
174 /* compute equivalent grey */
175 v = (brt - brtmin)/(brtmax - brtmin);
176 for (i = 0; i < 3; i++)
177 cgry[i] = v*upper[i] + (1.-v)*lower[i];
178 vv = 1.; /* check each limit */
179 for (i = 0; i < 3; i++)
180 if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
181 v = (lower[i]+CEPS - cgry[i])/(col[i] - cgry[i]);
182 if (v < vv) vv = v;
183 rflags |= CGAMUT_LOWER;
184 } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
185 v = (upper[i]-CEPS - cgry[i])/(col[i] - cgry[i]);
186 if (v < vv) vv = v;
187 rflags |= CGAMUT_UPPER;
188 }
189 if (rflags) /* desaturate to cube face */
190 for (i = 0; i < 3; i++)
191 col[i] = vv*col[i] + (1.-vv)*cgry[i];
192 return(rflags);
193 }
194
195
196 void
197 colortrans(c2, mat, c1) /* convert c1 by mat and put into c2 */
198 register COLOR c2;
199 register COLORMAT mat;
200 register COLOR c1;
201 {
202 COLOR cout;
203
204 cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
205 cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
206 cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];
207
208 copycolor(c2, cout);
209 }
210
211
212 void
213 multcolormat(m3, m2, m1) /* multiply m1 by m2 and put into m3 */
214 COLORMAT m3; /* m3 can be either m1 or m2 w/o harm */
215 COLORMAT m2, m1;
216 {
217 COLORMAT mt;
218 register int i, j;
219
220 for (i = 0; i < 3; i++)
221 for (j = 0; j < 3; j++)
222 mt[i][j] = m1[i][0]*m2[0][j] +
223 m1[i][1]*m2[1][j] +
224 m1[i][2]*m2[2][j] ;
225 cpcolormat(m3, mt);
226 }
227
228
229 void
230 compxyz2rgbmat(mat, pr) /* compute conversion from CIE to RGB space */
231 COLORMAT mat;
232 register RGBPRIMS pr;
233 {
234 double C_rD, C_gD, C_bD;
235
236 if (pr == stdprims) { /* can use xyz2rgbmat */
237 cpcolormat(mat, xyz2rgbmat);
238 return;
239 }
240 C_rD = (1./pr[WHT][CIEY]) *
241 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
242 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
243 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
244 C_gD = (1./pr[WHT][CIEY]) *
245 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
246 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
247 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
248 C_bD = (1./pr[WHT][CIEY]) *
249 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
250 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
251 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
252
253 mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
254 pr[BLU][CIEX]*pr[GRN][CIEY] +
255 pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
256 mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
257 pr[BLU][CIEX]*pr[GRN][CIEY] +
258 pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
259 mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
260 pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
261 mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
262 pr[BLU][CIEY]*pr[RED][CIEX] +
263 pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
264 mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
265 pr[RED][CIEX]*pr[BLU][CIEY] +
266 pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
267 mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
268 pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
269 mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
270 pr[RED][CIEY]*pr[GRN][CIEX] +
271 pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
272 mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
273 pr[GRN][CIEX]*pr[RED][CIEY] +
274 pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
275 mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
276 pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
277 }
278
279
280 void
281 comprgb2xyzmat(mat, pr) /* compute conversion from RGB to CIE space */
282 COLORMAT mat;
283 register RGBPRIMS pr;
284 {
285 double C_rD, C_gD, C_bD, D;
286
287 if (pr == stdprims) { /* can use rgb2xyzmat */
288 cpcolormat(mat, rgb2xyzmat);
289 return;
290 }
291 C_rD = (1./pr[WHT][CIEY]) *
292 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
293 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
294 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
295 C_gD = (1./pr[WHT][CIEY]) *
296 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
297 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
298 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
299 C_bD = (1./pr[WHT][CIEY]) *
300 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
301 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
302 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
303 D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
304 pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
305 pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
306 mat[0][0] = pr[RED][CIEX]*C_rD/D;
307 mat[0][1] = pr[GRN][CIEX]*C_gD/D;
308 mat[0][2] = pr[BLU][CIEX]*C_bD/D;
309 mat[1][0] = pr[RED][CIEY]*C_rD/D;
310 mat[1][1] = pr[GRN][CIEY]*C_gD/D;
311 mat[1][2] = pr[BLU][CIEY]*C_bD/D;
312 mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
313 mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
314 mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
315 }
316
317
318 void
319 comprgb2rgbmat(mat, pr1, pr2) /* compute conversion from RGB1 to RGB2 */
320 COLORMAT mat;
321 RGBPRIMS pr1, pr2;
322 {
323 COLORMAT pr1toxyz, xyztopr2;
324
325 if (pr1 == pr2) {
326 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
327 mat[0][1] = mat[0][2] = mat[1][0] =
328 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
329 return;
330 }
331 comprgb2xyzmat(pr1toxyz, pr1);
332 compxyz2rgbmat(xyztopr2, pr2);
333 /* combine transforms */
334 multcolormat(mat, pr1toxyz, xyztopr2);
335 }
336
337
338 void
339 compxyzWBmat(mat, wht1, wht2) /* CIE von Kries transform from wht1 to wht2 */
340 COLORMAT mat;
341 float wht1[2], wht2[2];
342 {
343 COLOR cw1, cw2;
344 if (XYEQ(wht1,wht2)) {
345 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
346 mat[0][1] = mat[0][2] = mat[1][0] =
347 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
348 return;
349 }
350 cw1[RED] = wht1[CIEX]/wht1[CIEY];
351 cw1[GRN] = 1.;
352 cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
353 colortrans(cw1, vkmat, cw1);
354 cw2[RED] = wht2[CIEX]/wht2[CIEY];
355 cw2[GRN] = 1.;
356 cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
357 colortrans(cw2, vkmat, cw2);
358 mat[0][0] = cw2[RED]/cw1[RED];
359 mat[1][1] = cw2[GRN]/cw1[GRN];
360 mat[2][2] = cw2[BLU]/cw1[BLU];
361 mat[0][1] = mat[0][2] = mat[1][0] =
362 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
363 multcolormat(mat, vkmat, mat);
364 multcolormat(mat, mat, ivkmat);
365 }
366
367
368 void
369 compxyz2rgbWBmat(mat, pr) /* von Kries conversion from CIE to RGB space */
370 COLORMAT mat;
371 RGBPRIMS pr;
372 {
373 COLORMAT wbmat;
374
375 compxyz2rgbmat(mat, pr);
376 if (XYEQ(pr[WHT],xyneu))
377 return;
378 compxyzWBmat(wbmat, xyneu, pr[WHT]);
379 multcolormat(mat, wbmat, mat);
380 }
381
382 void
383 comprgb2xyzWBmat(mat, pr) /* von Kries conversion from RGB to CIE space */
384 COLORMAT mat;
385 RGBPRIMS pr;
386 {
387 COLORMAT wbmat;
388
389 comprgb2xyzmat(mat, pr);
390 if (XYEQ(pr[WHT],xyneu))
391 return;
392 compxyzWBmat(wbmat, pr[WHT], xyneu);
393 multcolormat(mat, mat, wbmat);
394 }
395
396 void
397 comprgb2rgbWBmat(mat, pr1, pr2) /* von Kries conversion from RGB1 to RGB2 */
398 COLORMAT mat;
399 RGBPRIMS pr1, pr2;
400 {
401 COLORMAT pr1toxyz, xyztopr2, wbmat;
402
403 if (pr1 == pr2) {
404 mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
405 mat[0][1] = mat[0][2] = mat[1][0] =
406 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
407 return;
408 }
409 comprgb2xyzmat(pr1toxyz, pr1);
410 compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]);
411 compxyz2rgbmat(xyztopr2, pr2);
412 /* combine transforms */
413 multcolormat(mat, pr1toxyz, wbmat);
414 multcolormat(mat, mat, xyztopr2);
415 }