ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.26
Committed: Wed Sep 16 00:42:36 2020 UTC (4 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4
Changes since 2.25: +15 -9 lines
Log Message:
feat: Added functional check for XYZ color space to coloprimsOK()

File Contents

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