Back to index

radiance  4R0+20100331
spec_rgb.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: spec_rgb.c,v 2.20 2010/01/04 19:40:14 greg Exp $";
00003 #endif
00004 /*
00005  * Convert colors and spectral ranges.
00006  * Added von Kries white-balance calculations 10/01 (GW).
00007  *
00008  * Externals declared in color.h
00009  */
00010 
00011 #include "copyright.h"
00012 
00013 #include <stdio.h>
00014 #include <string.h>
00015 #include "color.h"
00016 
00017 #define CEPS  1e-4                 /* color epsilon */
00018 
00019 #define CEQ(v1,v2)   ((v1) <= (v2)+CEPS && (v2) <= (v1)+CEPS)
00020 
00021 #define XYEQ(c1,c2)  (CEQ((c1)[CIEX],(c2)[CIEX]) && CEQ((c1)[CIEY],(c2)[CIEY]))
00022 
00023 
00024 RGBPRIMS  stdprims = STDPRIMS;     /* standard primary chromaticities */
00025 
00026 COLOR  cblack = BLKCOLOR;          /* global black color */
00027 COLOR  cwhite = WHTCOLOR;          /* global white color */
00028 
00029 float  xyneu[2] = {1./3., 1./3.};  /* neutral xy chromaticities */
00030 
00031 /*
00032  *     The following table contains the CIE tristimulus integrals
00033  *  for X, Y, and Z.  The table is cumulative, so that
00034  *  each color coordinate integrates to 1.
00035  */
00036 
00037 #define  STARTWL     380           /* starting wavelength (nanometers) */
00038 #define  INCWL              10            /* wavelength increment */
00039 #define  NINC        40            /* # of values */
00040 
00041 static BYTE  chroma[3][NINC] = {
00042        {                                                /* X */
00043               0,   0,   0,   2,   6,   13,  22,  30,  36,  41,
00044               42,  43,  43,  44,  46,  52,  60,  71,  87,  106,
00045               128, 153, 178, 200, 219, 233, 243, 249, 252, 254,
00046               255, 255, 255, 255, 255, 255, 255, 255, 255, 255
00047        }, {                                             /* Y */
00048               0,   0,   0,   0,   0,   1,   2,   4,   7,   11,
00049               17,  24,  34,  48,  64,  84,  105, 127, 148, 169,
00050               188, 205, 220, 232, 240, 246, 250, 253, 254, 255,
00051               255, 255, 255, 255, 255, 255, 255, 255, 255, 255
00052        }, {                                             /* Z */
00053               0,   0,   2,   10,  32,  66,  118, 153, 191, 220,
00054               237, 246, 251, 253, 254, 255, 255, 255, 255, 255,
00055               255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
00056               255, 255, 255, 255, 255, 255, 255, 255, 255, 255
00057        }
00058 };
00059 
00060 COLORMAT  xyz2rgbmat = {           /* XYZ to RGB (no white balance) */
00061        {(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g)/CIE_C_rD,
00062         (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b)/CIE_C_rD,
00063         (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g)/CIE_C_rD},
00064        {(CIE_y_b - CIE_y_r - CIE_y_b*CIE_x_r + CIE_y_r*CIE_x_b)/CIE_C_gD,
00065         (CIE_x_r - CIE_x_b - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r)/CIE_C_gD,
00066         (CIE_x_b*CIE_y_r - CIE_x_r*CIE_y_b)/CIE_C_gD},
00067        {(CIE_y_r - CIE_y_g - CIE_y_r*CIE_x_g + CIE_y_g*CIE_x_r)/CIE_C_bD,
00068         (CIE_x_g - CIE_x_r - CIE_x_g*CIE_y_r + CIE_x_r*CIE_y_g)/CIE_C_bD,
00069         (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r)/CIE_C_bD}
00070 };
00071 
00072 COLORMAT  rgb2xyzmat = {           /* RGB to XYZ (no white balance) */
00073        {CIE_x_r*CIE_C_rD/CIE_D,CIE_x_g*CIE_C_gD/CIE_D,CIE_x_b*CIE_C_bD/CIE_D},
00074        {CIE_y_r*CIE_C_rD/CIE_D,CIE_y_g*CIE_C_gD/CIE_D,CIE_y_b*CIE_C_bD/CIE_D},
00075        {(1.-CIE_x_r-CIE_y_r)*CIE_C_rD/CIE_D,
00076         (1.-CIE_x_g-CIE_y_g)*CIE_C_gD/CIE_D,
00077         (1.-CIE_x_b-CIE_y_b)*CIE_C_bD/CIE_D}
00078 };
00079 
00080 COLORMAT  vkmat = {         /* Sharp primary matrix */
00081        { 1.2694, -0.0988, -0.1706},
00082        {-0.8364,  1.8006,  0.0357},
00083        { 0.0297, -0.0315,  1.0018}
00084 };
00085 
00086 COLORMAT  ivkmat = {        /* inverse Sharp primary matrix */
00087        { 0.8156,  0.0472,  0.1372},
00088        { 0.3791,  0.5769,  0.0440},
00089        {-0.0123,  0.0167,  0.9955}
00090 };
00091 
00092 
00093 void
00094 spec_rgb(                   /* compute RGB color from spectral range */
00095 COLOR  col,
00096 int  s,
00097 int  e
00098 )
00099 {
00100        COLOR  ciecolor;
00101 
00102        spec_cie(ciecolor, s, e);
00103        cie_rgb(col, ciecolor);
00104 }
00105 
00106 
00107 void
00108 spec_cie(                   /* compute a color from a spectral range */
00109 COLOR  col,          /* returned color */
00110 int  s,                     /* starting and ending wavelengths */
00111 int  e
00112 )
00113 {
00114        register int  i, d, r;
00115        
00116        s -= STARTWL;
00117        if (s < 0)
00118               s = 0;
00119 
00120        e -= STARTWL;
00121        if (e <= s) {
00122               col[CIEX] = col[CIEY] = col[CIEZ] = 0.0;
00123               return;
00124        }
00125        if (e >= INCWL*(NINC - 1))
00126               e = INCWL*(NINC - 1) - 1;
00127 
00128        d = e / INCWL;                     /* interpolate values */
00129        r = e % INCWL;
00130        for (i = 0; i < 3; i++)
00131               col[i] = chroma[i][d]*(INCWL - r) + chroma[i][d + 1]*r;
00132 
00133        d = s / INCWL;
00134        r = s % INCWL;
00135        for (i = 0; i < 3; i++)
00136               col[i] -= chroma[i][d]*(INCWL - r) - chroma[i][d + 1]*r;
00137 
00138        col[CIEX] = (col[CIEX] + 0.5) * (1./(256*INCWL));
00139        col[CIEY] = (col[CIEY] + 0.5) * (1./(256*INCWL));
00140        col[CIEZ] = (col[CIEZ] + 0.5) * (1./(256*INCWL));
00141 }
00142 
00143 
00144 void
00145 cie_rgb(                    /* convert CIE color to standard RGB */
00146 COLOR  rgb,
00147 COLOR  xyz
00148 )
00149 {
00150        colortrans(rgb, xyz2rgbmat, xyz);
00151        clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
00152 }
00153 
00154 
00155 int
00156 clipgamut(                  /* clip to gamut cube */
00157 COLOR  col,
00158 double  brt,
00159 int  gamut,
00160 COLOR  lower,
00161 COLOR  upper
00162 )
00163 {
00164        int  rflags = 0;
00165        double  brtmin, brtmax, v, vv;
00166        COLOR  cgry;
00167        register int  i;
00168                                    /* check for no check */
00169        if (gamut == 0) return(0);
00170                                    /* check brightness limits */
00171        brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
00172        if (gamut & CGAMUT_LOWER && brt < brtmin) {
00173               copycolor(col, lower);
00174               return(CGAMUT_LOWER);
00175        }
00176        brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
00177        if (gamut & CGAMUT_UPPER && brt > brtmax) {
00178               copycolor(col, upper);
00179               return(CGAMUT_UPPER);
00180        }
00181                                    /* compute equivalent grey */
00182        v = (brt - brtmin)/(brtmax - brtmin);
00183        for (i = 0; i < 3; i++)
00184               cgry[i] = v*upper[i] + (1.-v)*lower[i];
00185        vv = 1.;                    /* check each limit */
00186        for (i = 0; i < 3; i++)
00187               if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
00188                      v = (lower[i] - cgry[i])/(col[i] - cgry[i]);
00189                      if (v < vv) vv = v;
00190                      rflags |= CGAMUT_LOWER;
00191               } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
00192                      v = (upper[i] - cgry[i])/(col[i] - cgry[i]);
00193                      if (v < vv) vv = v;
00194                      rflags |= CGAMUT_UPPER;
00195               }
00196        if (rflags)                 /* desaturate to cube face */
00197               for (i = 0; i < 3; i++)
00198                      col[i] = vv*col[i] + (1.-vv)*cgry[i];
00199        return(rflags);
00200 }
00201 
00202 
00203 void
00204 colortrans(                 /* convert c1 by mat and put into c2 */
00205 register COLOR  c2,
00206 register COLORMAT  mat,
00207 register COLOR  c1
00208 )
00209 {
00210        COLOR  cout;
00211 
00212        cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
00213        cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
00214        cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];
00215 
00216        copycolor(c2, cout);
00217 }
00218 
00219 
00220 void
00221 multcolormat(               /* multiply m1 by m2 and put into m3 */
00222 COLORMAT  m3,               /* m3 can be either m1 or m2 w/o harm */
00223 COLORMAT  m2,
00224 COLORMAT  m1
00225 )
00226 {
00227        COLORMAT  mt;
00228        register int  i, j;
00229 
00230        for (i = 0; i < 3; i++)
00231               for (j = 0; j < 3; j++)
00232                      mt[i][j] =    m1[i][0]*m2[0][j] +
00233                                    m1[i][1]*m2[1][j] +
00234                                    m1[i][2]*m2[2][j] ;
00235        cpcolormat(m3, mt);
00236 }
00237 
00238 
00239 int
00240 colorprimsOK(               /* are color primaries reasonable? */
00241 RGBPRIMS  pr
00242 )
00243 {
00244        int    i, j;
00245        
00246        for (i = 0; i < 4; i++) {
00247               if ((pr[i][CIEX] <= -.5) | (pr[i][CIEY] <= -.5))
00248                      return(0);
00249               if ((pr[i][CIEX] >= 1.5) | (pr[i][CIEY] >= 1.5))
00250                      return(0);
00251               if (pr[i][CIEX] + pr[i][CIEY] >= 1.5)
00252                      return(0);
00253        }
00254        for (i = 0; i < 4; i++)
00255               for (j = i+1; j < 4; j++)
00256                      if (CEQ(pr[i][CIEX],pr[j][CIEX]) &&
00257                                    CEQ(pr[i][CIEY],pr[j][CIEY]))
00258                             return(0);
00259        return(1);
00260 }
00261 
00262 
00263 
00264 int
00265 compxyz2rgbmat(                    /* compute conversion from CIE to RGB space */
00266 COLORMAT  mat,
00267 register RGBPRIMS  pr
00268 )
00269 {
00270        double  C_rD, C_gD, C_bD;
00271 
00272        if (pr == stdprims) {       /* can use xyz2rgbmat */
00273               cpcolormat(mat, xyz2rgbmat);
00274               return(1);
00275        }
00276        if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
00277               return(0);
00278        C_rD = (1./pr[WHT][CIEY]) *
00279                      ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
00280                        pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
00281                 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
00282        if (CEQ(C_rD,0.))
00283               return(0);
00284        C_gD = (1./pr[WHT][CIEY]) *
00285                      ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
00286                        pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
00287                 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
00288        if (CEQ(C_gD,0.))
00289               return(0);
00290        C_bD = (1./pr[WHT][CIEY]) *
00291                      ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
00292                        pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
00293                 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
00294        if (CEQ(C_bD,0.))
00295               return(0);
00296        mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
00297                      pr[BLU][CIEX]*pr[GRN][CIEY] +
00298                      pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
00299        mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
00300                      pr[BLU][CIEX]*pr[GRN][CIEY] +
00301                      pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
00302        mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
00303                      pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
00304        mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
00305                      pr[BLU][CIEY]*pr[RED][CIEX] +
00306                      pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
00307        mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
00308                      pr[RED][CIEX]*pr[BLU][CIEY] +
00309                      pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
00310        mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
00311                      pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
00312        mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
00313                      pr[RED][CIEY]*pr[GRN][CIEX] +
00314                      pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
00315        mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
00316                      pr[GRN][CIEX]*pr[RED][CIEY] +
00317                      pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
00318        mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
00319                      pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
00320        return(1);
00321 }
00322 
00323 
00324 int
00325 comprgb2xyzmat(                    /* compute conversion from RGB to CIE space */
00326 COLORMAT  mat,
00327 register RGBPRIMS  pr
00328 )
00329 {
00330        double  C_rD, C_gD, C_bD, D;
00331 
00332        if (pr == stdprims) {       /* can use rgb2xyzmat */
00333               cpcolormat(mat, rgb2xyzmat);
00334               return(1);
00335        }
00336        if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
00337               return(0);
00338        C_rD = (1./pr[WHT][CIEY]) *
00339                      ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
00340                        pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
00341                 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
00342        C_gD = (1./pr[WHT][CIEY]) *
00343                      ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
00344                        pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
00345                 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
00346        C_bD = (1./pr[WHT][CIEY]) *
00347                      ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
00348                        pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
00349                 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
00350        D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
00351                      pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
00352                      pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
00353        if (CEQ(D,0.))
00354               return(0);
00355        mat[0][0] = pr[RED][CIEX]*C_rD/D;
00356        mat[0][1] = pr[GRN][CIEX]*C_gD/D;
00357        mat[0][2] = pr[BLU][CIEX]*C_bD/D;
00358        mat[1][0] = pr[RED][CIEY]*C_rD/D;
00359        mat[1][1] = pr[GRN][CIEY]*C_gD/D;
00360        mat[1][2] = pr[BLU][CIEY]*C_bD/D;
00361        mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
00362        mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
00363        mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
00364        return(1);
00365 }
00366 
00367 
00368 int
00369 comprgb2rgbmat(                    /* compute conversion from RGB1 to RGB2 */
00370 COLORMAT  mat,
00371 RGBPRIMS  pr1,
00372 RGBPRIMS  pr2
00373 )
00374 {
00375        COLORMAT  pr1toxyz, xyztopr2;
00376 
00377        if (pr1 == pr2) {
00378               mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
00379               mat[0][1] = mat[0][2] = mat[1][0] =
00380               mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
00381               return(1);
00382        }
00383        if (!comprgb2xyzmat(pr1toxyz, pr1))
00384               return(0);
00385        if (!compxyz2rgbmat(xyztopr2, pr2))
00386               return(0);
00387                             /* combine transforms */
00388        multcolormat(mat, pr1toxyz, xyztopr2);
00389        return(1);
00390 }
00391 
00392 
00393 int
00394 compxyzWBmat(               /* CIE von Kries transform from wht1 to wht2 */
00395 COLORMAT  mat,
00396 float  wht1[2],
00397 float  wht2[2]
00398 )
00399 {
00400        COLOR  cw1, cw2;
00401        if (XYEQ(wht1,wht2)) {
00402               mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
00403               mat[0][1] = mat[0][2] = mat[1][0] =
00404               mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
00405               return(1);
00406        }
00407        if (CEQ(wht1[CIEX],0.) | CEQ(wht1[CIEY],0.))
00408               return(0);
00409        cw1[RED] = wht1[CIEX]/wht1[CIEY];
00410        cw1[GRN] = 1.;
00411        cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
00412        colortrans(cw1, vkmat, cw1);
00413        if (CEQ(wht2[CIEX],0.) | CEQ(wht2[CIEY],0.))
00414               return(0);
00415        cw2[RED] = wht2[CIEX]/wht2[CIEY];
00416        cw2[GRN] = 1.;
00417        cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
00418        colortrans(cw2, vkmat, cw2);
00419        if (CEQ(cw1[RED],0.) | CEQ(cw1[GRN],0.) | CEQ(cw1[BLU],0.))
00420               return(0);
00421        mat[0][0] = cw2[RED]/cw1[RED];
00422        mat[1][1] = cw2[GRN]/cw1[GRN];
00423        mat[2][2] = cw2[BLU]/cw1[BLU];
00424        mat[0][1] = mat[0][2] = mat[1][0] =
00425        mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
00426        multcolormat(mat, vkmat, mat);
00427        multcolormat(mat, mat, ivkmat);
00428        return(1);
00429 }
00430 
00431 
00432 int
00433 compxyz2rgbWBmat(           /* von Kries conversion from CIE to RGB space */
00434 COLORMAT  mat,
00435 RGBPRIMS  pr
00436 )
00437 {
00438        COLORMAT      wbmat;
00439 
00440        if (!compxyz2rgbmat(mat, pr))
00441               return(0);
00442        if (XYEQ(pr[WHT],xyneu))
00443               return(1);
00444        if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
00445               return(0);
00446        multcolormat(mat, wbmat, mat);
00447        return(1);
00448 }
00449 
00450 int
00451 comprgb2xyzWBmat(           /* von Kries conversion from RGB to CIE space */
00452 COLORMAT  mat,
00453 RGBPRIMS  pr
00454 )
00455 {
00456        COLORMAT      wbmat;
00457        
00458        if (!comprgb2xyzmat(mat, pr))
00459               return(0);
00460        if (XYEQ(pr[WHT],xyneu))
00461               return(1);
00462        if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
00463               return(0);
00464        multcolormat(mat, mat, wbmat);
00465        return(1);
00466 }
00467 
00468 int
00469 comprgb2rgbWBmat(           /* von Kries conversion from RGB1 to RGB2 */
00470 COLORMAT  mat,
00471 RGBPRIMS  pr1,
00472 RGBPRIMS  pr2
00473 )
00474 {
00475        COLORMAT  pr1toxyz, xyztopr2, wbmat;
00476 
00477        if (pr1 == pr2) {
00478               mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
00479               mat[0][1] = mat[0][2] = mat[1][0] =
00480               mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
00481               return(1);
00482        }
00483        if (!comprgb2xyzmat(pr1toxyz, pr1))
00484               return(0);
00485        if (!compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]))
00486               return(0);
00487        if (!compxyz2rgbmat(xyztopr2, pr2))
00488               return(0);
00489                             /* combine transforms */
00490        multcolormat(mat, pr1toxyz, wbmat);
00491        multcolormat(mat, mat, xyztopr2);
00492        return(1);
00493 }