ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
(Generate patch)

Comparing ray/src/common/spec_rgb.c (file contents):
Revision 1.1 by greg, Sat Sep 22 10:44:56 1990 UTC vs.
Revision 2.16 by greg, Mon Oct 8 18:07:55 2007 UTC

# Line 1 | Line 1
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
4
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
# Line 37 | Line 57 | static BYTE  chroma[3][NINC] = {
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 < spec_rgb(col, s, e)             /* compute RGB color from spectral range */
73 < COLOR  col;
74 < int  s, e;
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  
# Line 49 | Line 104 | int  s, e;
104   }
105  
106  
107 < spec_cie(col, s, e)             /* compute a color from a spectral range */
108 < COLOR  col;             /* returned color */
109 < int  s, e;              /* starting and ending wavelengths */
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          
# Line 60 | Line 118 | int  s, e;             /* starting and ending wavelengths */
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  
# Line 73 | Line 135 | int  s, e;             /* starting and ending wavelengths */
135          for (i = 0; i < 3; i++)
136                  col[i] -= chroma[i][d]*(INCWL - r) - chroma[i][d + 1]*r;
137  
138 <        col[RED] = (col[RED] + 0.5) / (256*INCWL);
139 <        col[GRN] = (col[GRN] + 0.5) / (256*INCWL);
140 <        col[BLU] = (col[BLU] + 0.5) / (256*INCWL);
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 < cie_rgb(rgbcolor, ciecolor)             /* convert CIE to RGB (NTSC) */
145 < register COLOR  rgbcolor, ciecolor;
144 > void
145 > cie_rgb(                        /* convert CIE color to standard RGB */
146 > COLOR   rgb,
147 > COLOR  xyz
148 > )
149   {
150 <        static float  cmat[3][3] = {
151 <                1.73, -.48, -.26,
152 <                -.81, 1.65, -.02,
153 <                 .08, -.17, 1.28,
154 <        };
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 <        for (i = 0; i < 3; i++) {
203 <                rgbcolor[i] =   cmat[i][0]*ciecolor[0] +
204 <                                cmat[i][1]*ciecolor[1] +
205 <                                cmat[i][2]*ciecolor[2] ;
206 <                if (rgbcolor[i] < 0.0)
207 <                        rgbcolor[i] = 0.0;
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          }
252 +        return(1);
253 + }
254 +
255 +
256 +
257 + int
258 + compxyz2rgbmat(                 /* compute conversion from CIE to RGB space */
259 + COLORMAT  mat,
260 + register RGBPRIMS  pr
261 + )
262 + {
263 +        double  C_rD, C_gD, C_bD;
264 +
265 +        if (pr == stdprims) {   /* can use xyz2rgbmat */
266 +                cpcolormat(mat, xyz2rgbmat);
267 +                return(1);
268 +        }
269 +        if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
270 +                return(0);
271 +        C_rD = (1./pr[WHT][CIEY]) *
272 +                        ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
273 +                          pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
274 +                  pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
275 +        if (CEQ(C_rD,0.))
276 +                return(0);
277 +        C_gD = (1./pr[WHT][CIEY]) *
278 +                        ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
279 +                          pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
280 +                  pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
281 +        if (CEQ(C_gD,0.))
282 +                return(0);
283 +        C_bD = (1./pr[WHT][CIEY]) *
284 +                        ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
285 +                          pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
286 +                  pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
287 +        if (CEQ(C_bD,0.))
288 +                return(0);
289 +        mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
290 +                        pr[BLU][CIEX]*pr[GRN][CIEY] +
291 +                        pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
292 +        mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
293 +                        pr[BLU][CIEX]*pr[GRN][CIEY] +
294 +                        pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
295 +        mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
296 +                        pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
297 +        mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
298 +                        pr[BLU][CIEY]*pr[RED][CIEX] +
299 +                        pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
300 +        mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
301 +                        pr[RED][CIEX]*pr[BLU][CIEY] +
302 +                        pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
303 +        mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
304 +                        pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
305 +        mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
306 +                        pr[RED][CIEY]*pr[GRN][CIEX] +
307 +                        pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
308 +        mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
309 +                        pr[GRN][CIEX]*pr[RED][CIEY] +
310 +                        pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
311 +        mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
312 +                        pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
313 +        return(1);
314 + }
315 +
316 +
317 + int
318 + comprgb2xyzmat(                 /* compute conversion from RGB to CIE space */
319 + COLORMAT  mat,
320 + register RGBPRIMS  pr
321 + )
322 + {
323 +        double  C_rD, C_gD, C_bD, D;
324 +
325 +        if (pr == stdprims) {   /* can use rgb2xyzmat */
326 +                cpcolormat(mat, rgb2xyzmat);
327 +                return(1);
328 +        }
329 +        if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
330 +                return(0);
331 +        C_rD = (1./pr[WHT][CIEY]) *
332 +                        ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
333 +                          pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
334 +                  pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
335 +        C_gD = (1./pr[WHT][CIEY]) *
336 +                        ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
337 +                          pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
338 +                  pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
339 +        C_bD = (1./pr[WHT][CIEY]) *
340 +                        ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
341 +                          pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
342 +                  pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
343 +        D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
344 +                        pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
345 +                        pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
346 +        if (CEQ(D,0.))
347 +                return(0);
348 +        mat[0][0] = pr[RED][CIEX]*C_rD/D;
349 +        mat[0][1] = pr[GRN][CIEX]*C_gD/D;
350 +        mat[0][2] = pr[BLU][CIEX]*C_bD/D;
351 +        mat[1][0] = pr[RED][CIEY]*C_rD/D;
352 +        mat[1][1] = pr[GRN][CIEY]*C_gD/D;
353 +        mat[1][2] = pr[BLU][CIEY]*C_bD/D;
354 +        mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
355 +        mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
356 +        mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
357 +        return(1);
358 + }
359 +
360 +
361 + int
362 + comprgb2rgbmat(                 /* compute conversion from RGB1 to RGB2 */
363 + COLORMAT  mat,
364 + RGBPRIMS  pr1,
365 + RGBPRIMS  pr2
366 + )
367 + {
368 +        COLORMAT  pr1toxyz, xyztopr2;
369 +
370 +        if (pr1 == pr2) {
371 +                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
372 +                mat[0][1] = mat[0][2] = mat[1][0] =
373 +                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
374 +                return(1);
375 +        }
376 +        if (!comprgb2xyzmat(pr1toxyz, pr1))
377 +                return(0);
378 +        if (!compxyz2rgbmat(xyztopr2, pr2))
379 +                return(0);
380 +                                /* combine transforms */
381 +        multcolormat(mat, pr1toxyz, xyztopr2);
382 +        return(1);
383 + }
384 +
385 +
386 + int
387 + compxyzWBmat(                   /* CIE von Kries transform from wht1 to wht2 */
388 + COLORMAT  mat,
389 + float  wht1[2],
390 + float  wht2[2]
391 + )
392 + {
393 +        COLOR   cw1, cw2;
394 +        if (XYEQ(wht1,wht2)) {
395 +                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
396 +                mat[0][1] = mat[0][2] = mat[1][0] =
397 +                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
398 +                return(1);
399 +        }
400 +        if (CEQ(wht1[CIEX],0.) | CEQ(wht1[CIEY],0.))
401 +                return(0);
402 +        cw1[RED] = wht1[CIEX]/wht1[CIEY];
403 +        cw1[GRN] = 1.;
404 +        cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
405 +        colortrans(cw1, vkmat, cw1);
406 +        if (CEQ(wht2[CIEX],0.) | CEQ(wht2[CIEY],0.))
407 +                return(0);
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 +        if (CEQ(cw1[RED],0.) | CEQ(cw1[GRN],0.) | CEQ(cw1[BLU],0.))
413 +                return(0);
414 +        mat[0][0] = cw2[RED]/cw1[RED];
415 +        mat[1][1] = cw2[GRN]/cw1[GRN];
416 +        mat[2][2] = cw2[BLU]/cw1[BLU];
417 +        mat[0][1] = mat[0][2] = mat[1][0] =
418 +        mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
419 +        multcolormat(mat, vkmat, mat);
420 +        multcolormat(mat, mat, ivkmat);
421 +        return(1);
422 + }
423 +
424 +
425 + int
426 + compxyz2rgbWBmat(               /* von Kries conversion from CIE to RGB space */
427 + COLORMAT  mat,
428 + RGBPRIMS  pr
429 + )
430 + {
431 +        COLORMAT        wbmat;
432 +
433 +        if (!compxyz2rgbmat(mat, pr))
434 +                return(0);
435 +        if (XYEQ(pr[WHT],xyneu))
436 +                return(1);
437 +        if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
438 +                return(0);
439 +        multcolormat(mat, wbmat, mat);
440 +        return(1);
441 + }
442 +
443 + int
444 + comprgb2xyzWBmat(               /* von Kries conversion from RGB to CIE space */
445 + COLORMAT  mat,
446 + RGBPRIMS  pr
447 + )
448 + {
449 +        COLORMAT        wbmat;
450 +        
451 +        if (!comprgb2xyzmat(mat, pr))
452 +                return(0);
453 +        if (XYEQ(pr[WHT],xyneu))
454 +                return(1);
455 +        if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
456 +                return(0);
457 +        multcolormat(mat, mat, wbmat);
458 +        return(1);
459 + }
460 +
461 + int
462 + comprgb2rgbWBmat(               /* von Kries conversion from RGB1 to RGB2 */
463 + COLORMAT  mat,
464 + RGBPRIMS  pr1,
465 + RGBPRIMS  pr2
466 + )
467 + {
468 +        COLORMAT  pr1toxyz, xyztopr2, wbmat;
469 +
470 +        if (pr1 == pr2) {
471 +                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
472 +                mat[0][1] = mat[0][2] = mat[1][0] =
473 +                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
474 +                return(1);
475 +        }
476 +        if (!comprgb2xyzmat(pr1toxyz, pr1))
477 +                return(0);
478 +        if (!compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]))
479 +                return(0);
480 +        if (!compxyz2rgbmat(xyztopr2, pr2))
481 +                return(0);
482 +                                /* combine transforms */
483 +        multcolormat(mat, pr1toxyz, wbmat);
484 +        multcolormat(mat, mat, xyztopr2);
485 +        return(1);
486   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines