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.3 by greg, Fri Dec 14 12:48:32 1990 UTC vs.
Revision 2.12 by greg, Tue Feb 25 02:47:22 2003 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1990 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
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 "color.h"
14 + #include <string.h>
15  
16 + #define CEPS    1e-4                    /* color epsilon */
17 +
18 + #define CEQ(v1,v2)      ((v1) <= (v2)+CEPS && (v2) <= (v1)+CEPS)
19 +
20 + #define XYEQ(c1,c2)     (CEQ((c1)[CIEX],(c2)[CIEX]) && CEQ((c1)[CIEY],(c2)[CIEY]))
21 +
22 +
23 + RGBPRIMS  stdprims = STDPRIMS;  /* standard primary chromaticities */
24 +
25 + COLOR  cblack = BLKCOLOR;               /* global black color */
26 + COLOR  cwhite = WHTCOLOR;               /* global white color */
27 +
28 + float  xyneu[2] = {1./3., 1./3.};       /* neutral xy chromaticities */
29 +
30   /*
31   *      The following table contains the CIE tristimulus integrals
32   *  for X, Y, and Z.  The table is cumulative, so that
# Line 39 | Line 56 | static BYTE  chroma[3][NINC] = {
56          }
57   };
58  
59 < #ifdef  NTSC
60 < static float  xyz2rgbmat[3][3] = {      /* XYZ to RGB (NTSC) */
61 <        1.73, -.48, -.26,
62 <        -.81, 1.65, -.02,
63 <         .08, -.17, 1.28,
59 > COLORMAT  xyz2rgbmat = {                /* XYZ to RGB (no white balance) */
60 >        {(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g)/CIE_C_rD,
61 >         (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b)/CIE_C_rD,
62 >         (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g)/CIE_C_rD},
63 >        {(CIE_y_b - CIE_y_r - CIE_y_b*CIE_x_r + CIE_y_r*CIE_x_b)/CIE_C_gD,
64 >         (CIE_x_r - CIE_x_b - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r)/CIE_C_gD,
65 >         (CIE_x_b*CIE_y_r - CIE_x_r*CIE_y_b)/CIE_C_gD},
66 >        {(CIE_y_r - CIE_y_g - CIE_y_r*CIE_x_g + CIE_y_g*CIE_x_r)/CIE_C_bD,
67 >         (CIE_x_g - CIE_x_r - CIE_x_g*CIE_y_r + CIE_x_r*CIE_y_g)/CIE_C_bD,
68 >         (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r)/CIE_C_bD}
69   };
70 < #else
71 < static float xyz2rgbmat[3][3] = {       /* XYZ to RGB (color monitor) */
72 <         2.739, -1.145, -.424,
73 <        -1.119,  2.029,  .033,
74 <          .138,  -.333, 1.105,
70 >
71 > COLORMAT  rgb2xyzmat = {                /* RGB to XYZ (no white balance) */
72 >        {CIE_x_r*CIE_C_rD/CIE_D,CIE_x_g*CIE_C_gD/CIE_D,CIE_x_b*CIE_C_bD/CIE_D},
73 >        {CIE_y_r*CIE_C_rD/CIE_D,CIE_y_g*CIE_C_gD/CIE_D,CIE_y_b*CIE_C_bD/CIE_D},
74 >        {(1.-CIE_x_r-CIE_y_r)*CIE_C_rD/CIE_D,
75 >         (1.-CIE_x_g-CIE_y_g)*CIE_C_gD/CIE_D,
76 >         (1.-CIE_x_b-CIE_y_b)*CIE_C_bD/CIE_D}
77   };
54 – #endif
78  
79 + COLORMAT  vkmat = {             /* Sharp primary matrix */
80 +        { 1.2694, -0.0988, -0.1706},
81 +        {-0.8364,  1.8006,  0.0357},
82 +        { 0.0297, -0.0315,  1.0018}
83 + };
84  
85 + COLORMAT  ivkmat = {            /* inverse Sharp primary matrix */
86 +        { 0.8156,  0.0472,  0.1372},
87 +        { 0.3791,  0.5769,  0.0440},
88 +        {-0.0123,  0.0167,  0.9955}
89 + };
90  
91 +
92 + void
93   spec_rgb(col, s, e)             /* compute RGB color from spectral range */
94   COLOR  col;
95   int  s, e;
# Line 66 | Line 101 | int  s, e;
101   }
102  
103  
104 + void
105   spec_cie(col, s, e)             /* compute a color from a spectral range */
106   COLOR  col;             /* returned color */
107   int  s, e;              /* starting and ending wavelengths */
# Line 78 | Line 114 | int  s, e;             /* starting and ending wavelengths */
114  
115          e -= STARTWL;
116          if (e <= s) {
117 <                col[RED] = col[GRN] = col[BLU] = 0.0;
117 >                col[CIEX] = col[CIEY] = col[CIEZ] = 0.0;
118                  return;
119          }
120          if (e >= INCWL*(NINC - 1))
# Line 94 | Line 130 | int  s, e;             /* starting and ending wavelengths */
130          for (i = 0; i < 3; i++)
131                  col[i] -= chroma[i][d]*(INCWL - r) - chroma[i][d + 1]*r;
132  
133 <        col[RED] = (col[RED] + 0.5) / (256*INCWL);
134 <        col[GRN] = (col[GRN] + 0.5) / (256*INCWL);
135 <        col[BLU] = (col[BLU] + 0.5) / (256*INCWL);
133 >        col[CIEX] = (col[CIEX] + 0.5) * (1./(256*INCWL));
134 >        col[CIEY] = (col[CIEY] + 0.5) * (1./(256*INCWL));
135 >        col[CIEZ] = (col[CIEZ] + 0.5) * (1./(256*INCWL));
136   }
137  
138  
139 < cie_rgb(rgbcolor, ciecolor)             /* convert CIE to RGB */
140 < register COLOR  rgbcolor, ciecolor;
139 > void
140 > cie_rgb(rgb, xyz)               /* convert CIE color to standard RGB */
141 > COLOR   rgb;
142 > COLOR  xyz;
143   {
144 +        colortrans(rgb, xyz2rgbmat, xyz);
145 +        clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
146 + }
147 +
148 +
149 + int
150 + clipgamut(col, brt, gamut, lower, upper)        /* clip to gamut cube */
151 + COLOR  col;
152 + double  brt;
153 + int  gamut;
154 + COLOR  lower, upper;
155 + {
156 +        int  rflags = 0;
157 +        double  brtmin, brtmax, v, vv;
158 +        COLOR  cgry;
159          register int  i;
160 +                                        /* check for no check */
161 +        if (gamut == 0) return(0);
162 +                                        /* check brightness limits */
163 +        brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
164 +        if (gamut & CGAMUT_LOWER && brt < brtmin) {
165 +                copycolor(col, lower);
166 +                return(CGAMUT_LOWER);
167 +        }
168 +        brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
169 +        if (gamut & CGAMUT_UPPER && brt > brtmax) {
170 +                copycolor(col, upper);
171 +                return(CGAMUT_UPPER);
172 +        }
173 +                                        /* compute equivalent grey */
174 +        v = (brt - brtmin)/(brtmax - brtmin);
175 +        for (i = 0; i < 3; i++)
176 +                cgry[i] = v*upper[i] + (1.-v)*lower[i];
177 +        vv = 1.;                        /* check each limit */
178 +        for (i = 0; i < 3; i++)
179 +                if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
180 +                        v = (lower[i]+CEPS - cgry[i])/(col[i] - cgry[i]);
181 +                        if (v < vv) vv = v;
182 +                        rflags |= CGAMUT_LOWER;
183 +                } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
184 +                        v = (upper[i]-CEPS - cgry[i])/(col[i] - cgry[i]);
185 +                        if (v < vv) vv = v;
186 +                        rflags |= CGAMUT_UPPER;
187 +                }
188 +        if (rflags)                     /* desaturate to cube face */
189 +                for (i = 0; i < 3; i++)
190 +                        col[i] = vv*col[i] + (1.-vv)*cgry[i];
191 +        return(rflags);
192 + }
193  
194 <        for (i = 0; i < 3; i++) {
195 <                rgbcolor[i] =   xyz2rgbmat[i][0]*ciecolor[0] +
196 <                                xyz2rgbmat[i][1]*ciecolor[1] +
197 <                                xyz2rgbmat[i][2]*ciecolor[2] ;
198 <                if (rgbcolor[i] < 0.0)
199 <                        rgbcolor[i] = 0.0;
194 >
195 > void
196 > colortrans(c2, mat, c1)         /* convert c1 by mat and put into c2 */
197 > register COLOR  c2;
198 > register COLORMAT  mat;
199 > register COLOR  c1;
200 > {
201 >        COLOR   cout;
202 >
203 >        cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
204 >        cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
205 >        cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];
206 >
207 >        copycolor(c2, cout);
208 > }
209 >
210 >
211 > void
212 > multcolormat(m3, m2, m1)        /* multiply m1 by m2 and put into m3 */
213 > COLORMAT  m3;                   /* m3 can be either m1 or m2 w/o harm */
214 > COLORMAT  m2, m1;
215 > {
216 >        COLORMAT  mt;
217 >        register int  i, j;
218 >
219 >        for (i = 0; i < 3; i++)
220 >                for (j = 0; j < 3; j++)
221 >                        mt[i][j] =      m1[i][0]*m2[0][j] +
222 >                                        m1[i][1]*m2[1][j] +
223 >                                        m1[i][2]*m2[2][j] ;
224 >        cpcolormat(m3, mt);
225 > }
226 >
227 >
228 > void
229 > compxyz2rgbmat(mat, pr)         /* compute conversion from CIE to RGB space */
230 > COLORMAT  mat;
231 > register RGBPRIMS  pr;
232 > {
233 >        double  C_rD, C_gD, C_bD;
234 >
235 >        if (pr == stdprims) {   /* can use xyz2rgbmat */
236 >                cpcolormat(mat, xyz2rgbmat);
237 >                return;
238          }
239 +        C_rD = (1./pr[WHT][CIEY]) *
240 +                        ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
241 +                          pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
242 +                  pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
243 +        C_gD = (1./pr[WHT][CIEY]) *
244 +                        ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
245 +                          pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
246 +                  pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
247 +        C_bD = (1./pr[WHT][CIEY]) *
248 +                        ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
249 +                          pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
250 +                  pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
251 +
252 +        mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
253 +                        pr[BLU][CIEX]*pr[GRN][CIEY] +
254 +                        pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
255 +        mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
256 +                        pr[BLU][CIEX]*pr[GRN][CIEY] +
257 +                        pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
258 +        mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
259 +                        pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
260 +        mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
261 +                        pr[BLU][CIEY]*pr[RED][CIEX] +
262 +                        pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
263 +        mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
264 +                        pr[RED][CIEX]*pr[BLU][CIEY] +
265 +                        pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
266 +        mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
267 +                        pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
268 +        mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
269 +                        pr[RED][CIEY]*pr[GRN][CIEX] +
270 +                        pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
271 +        mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
272 +                        pr[GRN][CIEX]*pr[RED][CIEY] +
273 +                        pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
274 +        mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
275 +                        pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
276 + }
277 +
278 +
279 + void
280 + comprgb2xyzmat(mat, pr)         /* compute conversion from RGB to CIE space */
281 + COLORMAT  mat;
282 + register RGBPRIMS  pr;
283 + {
284 +        double  C_rD, C_gD, C_bD, D;
285 +
286 +        if (pr == stdprims) {   /* can use rgb2xyzmat */
287 +                cpcolormat(mat, rgb2xyzmat);
288 +                return;
289 +        }
290 +        C_rD = (1./pr[WHT][CIEY]) *
291 +                        ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
292 +                          pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
293 +                  pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
294 +        C_gD = (1./pr[WHT][CIEY]) *
295 +                        ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
296 +                          pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
297 +                  pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
298 +        C_bD = (1./pr[WHT][CIEY]) *
299 +                        ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
300 +                          pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
301 +                  pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
302 +        D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
303 +                        pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
304 +                        pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
305 +        mat[0][0] = pr[RED][CIEX]*C_rD/D;
306 +        mat[0][1] = pr[GRN][CIEX]*C_gD/D;
307 +        mat[0][2] = pr[BLU][CIEX]*C_bD/D;
308 +        mat[1][0] = pr[RED][CIEY]*C_rD/D;
309 +        mat[1][1] = pr[GRN][CIEY]*C_gD/D;
310 +        mat[1][2] = pr[BLU][CIEY]*C_bD/D;
311 +        mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
312 +        mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
313 +        mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
314 + }
315 +
316 +
317 + void
318 + comprgb2rgbmat(mat, pr1, pr2)   /* compute conversion from RGB1 to RGB2 */
319 + COLORMAT  mat;
320 + RGBPRIMS  pr1, pr2;
321 + {
322 +        COLORMAT  pr1toxyz, xyztopr2;
323 +
324 +        if (pr1 == pr2) {
325 +                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
326 +                mat[0][1] = mat[0][2] = mat[1][0] =
327 +                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
328 +                return;
329 +        }
330 +        comprgb2xyzmat(pr1toxyz, pr1);
331 +        compxyz2rgbmat(xyztopr2, pr2);
332 +                                /* combine transforms */
333 +        multcolormat(mat, pr1toxyz, xyztopr2);
334 + }
335 +
336 +
337 + void
338 + compxyzWBmat(mat, wht1, wht2)   /* CIE von Kries transform from wht1 to wht2 */
339 + COLORMAT  mat;
340 + float  wht1[2], wht2[2];
341 + {
342 +        COLOR   cw1, cw2;
343 +        if (XYEQ(wht1,wht2)) {
344 +                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
345 +                mat[0][1] = mat[0][2] = mat[1][0] =
346 +                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
347 +                return;
348 +        }
349 +        cw1[RED] = wht1[CIEX]/wht1[CIEY];
350 +        cw1[GRN] = 1.;
351 +        cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
352 +        colortrans(cw1, vkmat, cw1);
353 +        cw2[RED] = wht2[CIEX]/wht2[CIEY];
354 +        cw2[GRN] = 1.;
355 +        cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
356 +        colortrans(cw2, vkmat, cw2);
357 +        mat[0][0] = cw2[RED]/cw1[RED];
358 +        mat[1][1] = cw2[GRN]/cw1[GRN];
359 +        mat[2][2] = cw2[BLU]/cw1[BLU];
360 +        mat[0][1] = mat[0][2] = mat[1][0] =
361 +        mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
362 +        multcolormat(mat, vkmat, mat);
363 +        multcolormat(mat, mat, ivkmat);
364 + }
365 +
366 +
367 + void
368 + compxyz2rgbWBmat(mat, pr)       /* von Kries conversion from CIE to RGB space */
369 + COLORMAT  mat;
370 + RGBPRIMS  pr;
371 + {
372 +        COLORMAT        wbmat;
373 +
374 +        compxyz2rgbmat(mat, pr);
375 +        if (XYEQ(pr[WHT],xyneu))
376 +                return;
377 +        compxyzWBmat(wbmat, xyneu, pr[WHT]);
378 +        multcolormat(mat, wbmat, mat);
379 + }
380 +
381 + void
382 + comprgb2xyzWBmat(mat, pr)       /* von Kries conversion from RGB to CIE space */
383 + COLORMAT  mat;
384 + RGBPRIMS  pr;
385 + {
386 +        COLORMAT        wbmat;
387 +        
388 +        comprgb2xyzmat(mat, pr);
389 +        if (XYEQ(pr[WHT],xyneu))
390 +                return;
391 +        compxyzWBmat(wbmat, pr[WHT], xyneu);
392 +        multcolormat(mat, mat, wbmat);
393 + }
394 +
395 + void
396 + comprgb2rgbWBmat(mat, pr1, pr2) /* von Kries conversion from RGB1 to RGB2 */
397 + COLORMAT  mat;
398 + RGBPRIMS  pr1, pr2;
399 + {
400 +        COLORMAT  pr1toxyz, xyztopr2, wbmat;
401 +
402 +        if (pr1 == pr2) {
403 +                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
404 +                mat[0][1] = mat[0][2] = mat[1][0] =
405 +                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
406 +                return;
407 +        }
408 +        comprgb2xyzmat(pr1toxyz, pr1);
409 +        compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]);
410 +        compxyz2rgbmat(xyztopr2, pr2);
411 +                                /* combine transforms */
412 +        multcolormat(mat, pr1toxyz, wbmat);
413 +        multcolormat(mat, mat, xyztopr2);
414   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines