home *** CD-ROM | disk | FTP | other *** search
/ Online Praxis 1996 April / OP4_96.ISO / surfen / astro / earthv / specrend.c < prev    next >
C/C++ Source or Header  |  1996-08-16  |  12KB  |  285 lines
  1. /*
  2.                      Colour Rendering of Spectra
  3.  
  4.                  by John Walker, kelvin@fourmilab.ch
  5.                        http://www.fourmilab.ch/
  6.  
  7.                 This program is in the public domain.
  8.  
  9.     For complete information about the techniques employed in this
  10.     program, see the World-Wide Web document:
  11.  
  12.              http://www.fourmilab.ch/documents/specrend/
  13.  
  14. */
  15.  
  16. #include <stdio.h>
  17. #include <math.h>
  18.  
  19. #define CLAMP(v,l,h)    ((v)<(l) ? (l) : (v) > (h) ? (h) : (v))
  20.  
  21. /* A colour system is defined by the CIE x and y coordinates of its
  22.    three primary illuminants and the x and y coordinates of the white
  23.    point. */
  24.  
  25. struct colourSystem {
  26.     double xRed, yRed,
  27.            xGreen, yGreen,
  28.            xBlue, yBlue,
  29.            xWhite, yWhite;
  30. };
  31.  
  32. /* White point chromaticities. */
  33.  
  34. #define IlluminantC     0.3101, 0.3162  /* For NTSC television */
  35. #define IlluminantD65   0.3127, 0.3291  /* For EBU and SMPTE */
  36.  
  37. static struct colourSystem
  38.                   /* xRed   yRed  xGreen yGreen  xBlue  yBlue  White point */
  39.     NTSCsystem  =  { 0.67,  0.33,  0.21,  0.71,  0.14,  0.08,  IlluminantC   },
  40.     EBUsystem   =  { 0.64,  0.33,  0.29,  0.60,  0.15,  0.06,  IlluminantD65 },
  41.     SMPTEsystem =  { 0.630, 0.340, 0.310, 0.595, 0.155, 0.070, IlluminantD65 };
  42.  
  43. /*                             XYZ_TO_RGB
  44.  
  45.     Given an additive tricolour system CS, defined by the  CIE  x  and  y
  46.     chromaticities  of  its  three  primaries (z is derived trivially as
  47.     1-(x+y)), and a desired chromaticity (XC,  YC,  ZC)  in  CIE  space,
  48.     determine  the  contribution of each primary in a linear combination
  49.     which  sums  to  the  desired  chromaticity.    If   the   requested
  50.     chromaticity falls outside the Maxwell triangle (colour gamut) formed
  51.     by the three primaries, one of the  r,  g,  or  b  weights  will  be
  52.     negative.   Use  inside_gamut()  to  test  for  a  valid  colour  and
  53.     constrain_rgb() to desaturate an outside-gamut colour to the  closest
  54.     representation within the available gamut. */
  55.  
  56.  
  57. void xyz_to_rgb(struct colourSystem *cs,
  58.                 double xc, double yc, double zc,
  59.                 double *r, double *g, double *b)
  60. {
  61.     double xr, yr, zr, xg, yg, zg, xb, yb, zb, d;
  62.  
  63.     xr = cs->xRed;    yr = cs->yRed;    zr = 1 - (xr + yr);
  64.     xg = cs->xGreen;  yg = cs->yGreen;  zg = 1 - (xg + yg);
  65.     xb = cs->xBlue;   yb = cs->yBlue;   zb = 1 - (xb + yb);
  66.     d = xr*yg*zb - xg*yr*zb - xr*yb*zg + xb*yr*zg + xg*yb*zr - xb*yg*zr;
  67.  
  68.     *r = (-xg*yc*zb + xc*yg*zb + xg*yb*zc - xb*yg*zc - xc*yb*zg + xb*yc*zg) / d;
  69.  
  70.     *g = (xr*yc*zb - xc*yr*zb - xr*yb*zc + xb*yr*zc + xc*yb*zr - xb*yc*zr) / d;
  71.  
  72.     *b = (xr*yg*zc - xg*yr*zc - xr*yc*zg + xc*yr*zg + xg*yc*zr - xc*yg*zr) / d;
  73. }
  74.  
  75. /*                            INSIDE_GAMUT
  76.  
  77.     Test  whether  a requested colour is within the gamut achievable with
  78.     the primaries of the current colour system.  This amounts  simply  to
  79.     testing whether all the primary weights are non-negative. */
  80.  
  81. int inside_gamut(double r, double g, double b)
  82. {
  83.     return r >= 0 && g >= 0 && b >= 0;
  84. }
  85.  
  86. /*                           CONSTRAIN_RGB
  87.  
  88.     If  the  requested RGB shade contains a negative weight for one of
  89.     the primaries, it lies outside the colour gamut accessible from the
  90.     given triple of primaries.  Desaturate it by mixing with the white
  91.     point of the colour system so as to reduce  the  primary  with  the
  92.     negative weight to  zero.   This  is  equivalent  to  finding  the
  93.     intersection  on the CIE diagram of a line drawn between the white
  94.     point and the  requested  colour  with  the  edge  of  the  Maxwell
  95.     triangle formed by the three primaries. */
  96.  
  97. int constrain_rgb(struct colourSystem *cs,
  98.                   double *x, double *y, double *z,
  99.                   double *r, double *g, double *b)
  100. {
  101.     /* Is the contribution of one of the primaries negative ? */
  102.  
  103.     if (!inside_gamut(*r, *g, *b)) {
  104.         double par, wr, wg, wb;
  105.  
  106.         /* Yes.  Find the RGB mixing weights of the white point (we
  107.            assume the white point is in the gamut!). */
  108.  
  109.         xyz_to_rgb(cs,
  110.                    cs->xWhite, cs->yWhite, 1 - (cs->xWhite + cs->yWhite),
  111.                    &wr, &wg, &wb);
  112.  
  113.         /* Find the primary with negative weight and calculate the
  114.            parameter of the point on the vector from the white point
  115.            to the original requested colour in RGB space. */
  116.  
  117.         if (*r < *g && *r < *b) {
  118.             par = wr / (wr - *r);
  119.         } else if (*g < *r && *g < *b) {
  120.             par = wg / (wg - *g);
  121.         } else {
  122.             par = wb / (wb - *b);
  123.         }
  124.  
  125.         /* Since XYZ space is a linear transformation of RGB space, we
  126.            can find the XYZ space coordinates of the point where the
  127.            edge of the gamut intersects the vector from the white point
  128.            to the original colour by multiplying the parameter in RGB
  129.            space by the difference vector in XYZ space. */
  130.  
  131.         *x = CLAMP(cs->xWhite + par * (*x - cs->xWhite), 0, 1);
  132.         *y = CLAMP(cs->yWhite + par * (*y - cs->yWhite), 0, 1);
  133.         *z = CLAMP(1 - (*x + *y), 0, 1);
  134.  
  135.         /* Now finally calculate the gamut-constrained RGB weights. */
  136.  
  137.         *r = CLAMP(wr + par * (*r - wr), 0, 1);
  138.         *g = CLAMP(wg + par * (*g - wg), 0, 1);
  139.         *b = CLAMP(wb + par * (*b - wb), 0, 1);
  140.         return 1;                     /* Colour modified to fit RGB gamut */
  141.     }
  142.     return 0;                         /* Colour within RGB gamut */
  143. }
  144.  
  145. /*                          SPECTRUM_TO_XYZ
  146.  
  147.     Calculate the CIE X, Y, and Z coordinates corresponding to a light
  148.     source  with  spectral  distribution   given   by   the   function
  149.     SPEC_INTENS,  which is called with a series of wavelengths between
  150.     380 and 780 nm  (the  argument  is  expressed  in  meters),  which
  151.     returns  emittance  at  that  wavelength  in arbitrary units.  The
  152.     chromaticity coordinates of the spectrum are returned in the x, y,
  153.     and z arguments which respect the identity:
  154.  
  155.             x + y + z = 1.
  156. */
  157.  
  158. void spectrum_to_xyz(double (*spec_intens)(double wavelength),
  159.                      double *x, double *y, double *z)
  160. {
  161.     int i;
  162.     double lambda, X = 0, Y = 0, Z = 0, XYZ;
  163.  
  164.     /* CIE colour matching functions xBar, yBar, and zBar for
  165.        wavelengths from 380 through 780 nanometers, every 5
  166.        nanometers.  For a wavelength lambda in this range:
  167.  
  168.             cie_colour_match[(lambda - 380) / 5][0] = xBar
  169.             cie_colour_match[(lambda - 380) / 5][1] = yBar
  170.             cie_colour_match[(lambda - 380) / 5][2] = zBar
  171.  
  172.        To  save  memory,  this  table can be declared as floats rather
  173.        than doubles; (IEEE)  float  has  enough  significant  bits  to
  174.        represent  the values.  It's declared as a double here to avoid
  175.        warnings about "conversion between floating-point  types"  from
  176.        certain persnickety compilers. */
  177.  
  178.     static double cie_colour_match[81][3] = {
  179.         {0.0014,0.0000,0.0065}, {0.0022,0.0001,0.0105}, {0.0042,0.0001,0.0201},
  180.         {0.0076,0.0002,0.0362}, {0.0143,0.0004,0.0679}, {0.0232,0.0006,0.1102},
  181.         {0.0435,0.0012,0.2074}, {0.0776,0.0022,0.3713}, {0.1344,0.0040,0.6456},
  182.         {0.2148,0.0073,1.0391}, {0.2839,0.0116,1.3856}, {0.3285,0.0168,1.6230},
  183.         {0.3483,0.0230,1.7471}, {0.3481,0.0298,1.7826}, {0.3362,0.0380,1.7721},
  184.         {0.3187,0.0480,1.7441}, {0.2908,0.0600,1.6692}, {0.2511,0.0739,1.5281},
  185.         {0.1954,0.0910,1.2876}, {0.1421,0.1126,1.0419}, {0.0956,0.1390,0.8130},
  186.         {0.0580,0.1693,0.6162}, {0.0320,0.2080,0.4652}, {0.0147,0.2586,0.3533},
  187.         {0.0049,0.3230,0.2720}, {0.0024,0.4073,0.2123}, {0.0093,0.5030,0.1582},
  188.         {0.0291,0.6082,0.1117}, {0.0633,0.7100,0.0782}, {0.1096,0.7932,0.0573},
  189.         {0.1655,0.8620,0.0422}, {0.2257,0.9149,0.0298}, {0.2904,0.9540,0.0203},
  190.         {0.3597,0.9803,0.0134}, {0.4334,0.9950,0.0087}, {0.5121,1.0000,0.0057},
  191.         {0.5945,0.9950,0.0039}, {0.6784,0.9786,0.0027}, {0.7621,0.9520,0.0021},
  192.         {0.8425,0.9154,0.0018}, {0.9163,0.8700,0.0017}, {0.9786,0.8163,0.0014},
  193.         {1.0263,0.7570,0.0011}, {1.0567,0.6949,0.0010}, {1.0622,0.6310,0.0008},
  194.         {1.0456,0.5668,0.0006}, {1.0026,0.5030,0.0003}, {0.9384,0.4412,0.0002},
  195.         {0.8544,0.3810,0.0002}, {0.7514,0.3210,0.0001}, {0.6424,0.2650,0.0000},
  196.         {0.5419,0.2170,0.0000}, {0.4479,0.1750,0.0000}, {0.3608,0.1382,0.0000},
  197.         {0.2835,0.1070,0.0000}, {0.2187,0.0816,0.0000}, {0.1649,0.0610,0.0000},
  198.         {0.1212,0.0446,0.0000}, {0.0874,0.0320,0.0000}, {0.0636,0.0232,0.0000},
  199.         {0.0468,0.0170,0.0000}, {0.0329,0.0119,0.0000}, {0.0227,0.0082,0.0000},
  200.         {0.0158,0.0057,0.0000}, {0.0114,0.0041,0.0000}, {0.0081,0.0029,0.0000},
  201.         {0.0058,0.0021,0.0000}, {0.0041,0.0015,0.0000}, {0.0029,0.0010,0.0000},
  202.         {0.0020,0.0007,0.0000}, {0.0014,0.0005,0.0000}, {0.0010,0.0004,0.0000},
  203.         {0.0007,0.0002,0.0000}, {0.0005,0.0002,0.0000}, {0.0003,0.0001,0.0000},
  204.         {0.0002,0.0001,0.0000}, {0.0002,0.0001,0.0000}, {0.0001,0.0000,0.0000},
  205.         {0.0001,0.0000,0.0000}, {0.0001,0.0000,0.0000}, {0.0000,0.0000,0.0000}
  206.     };
  207.  
  208.     for (i = 0, lambda = 380; lambda < 780.1; i++, lambda += 5) {
  209.         double Me;
  210.  
  211.         Me = (*spec_intens)(lambda);
  212.         X += Me * cie_colour_match[i][0];
  213.         Y += Me * cie_colour_match[i][1];
  214.         Z += Me * cie_colour_match[i][2];
  215.     }
  216.     XYZ = (X + Y + Z);
  217.     *x = X / XYZ;
  218.     *y = Y / XYZ;
  219.     *z = Z / XYZ;
  220. }
  221.  
  222. /*                            BB_SPECTRUM
  223.  
  224.     Calculate, by Planck's radiation law, the emittance of a black body
  225.     of temperature bbTemp at the given wavelength (in meters).  */
  226.  
  227. double bbTemp = 5000;                 /* Hidden temperature argument
  228.                                          to BB_SPECTRUM. */
  229. double bb_spectrum(double wavelength)
  230. {
  231.     double wlm = wavelength * 1e-9;   /* Wavelength in meters */
  232.  
  233.     return (3.74183e-16 * pow(wlm, -5.0)) /
  234.            (exp(1.4388e-2 / (wlm * bbTemp)) - 1.0);
  235. }
  236.  
  237. /*  Built-in test program which displays the x, y, and Z and RGB
  238.     values for black body spectra from 1000 to 10000 degrees kelvin.
  239.     When run, this program should produce the following output:
  240.  
  241.     Temperature       x      y      z       R     G     B
  242.     -----------    ------ ------ ------   ----- ----- -----
  243.        1000 K      0.6528 0.3444 0.0028   0.988 0.012 0.000 (Approximation)
  244.        1500 K      0.5857 0.3931 0.0212   0.808 0.192 0.000 (Approximation)
  245.        2000 K      0.5267 0.4133 0.0600   0.684 0.302 0.014
  246.        2500 K      0.4770 0.4137 0.1093   0.558 0.368 0.074
  247.        3000 K      0.4369 0.4041 0.1590   0.464 0.398 0.138
  248.        3500 K      0.4053 0.3907 0.2040   0.394 0.408 0.198
  249.        4000 K      0.3805 0.3768 0.2428   0.341 0.409 0.250
  250.        4500 K      0.3608 0.3636 0.2756   0.301 0.404 0.294
  251.        5000 K      0.3451 0.3516 0.3032   0.271 0.397 0.332
  252.        5500 K      0.3325 0.3411 0.3265   0.246 0.390 0.364
  253.        6000 K      0.3221 0.3318 0.3461   0.227 0.382 0.391
  254.        6500 K      0.3135 0.3237 0.3628   0.212 0.374 0.414
  255.        7000 K      0.3064 0.3166 0.3770   0.199 0.367 0.434
  256.        7500 K      0.3004 0.3103 0.3893   0.188 0.361 0.451
  257.        8000 K      0.2952 0.3048 0.4000   0.179 0.355 0.466
  258.        8500 K      0.2908 0.3000 0.4093   0.172 0.350 0.479
  259.        9000 K      0.2869 0.2956 0.4174   0.165 0.345 0.490
  260.        9500 K      0.2836 0.2918 0.4246   0.160 0.340 0.500
  261.       10000 K      0.2807 0.2884 0.4310   0.155 0.336 0.509
  262. */
  263.  
  264. int main()
  265. { 
  266.     double t, x, y, z, r, g, b;
  267.     struct colourSystem *cs = &SMPTEsystem;
  268.  
  269.     printf("Temperature       x      y      z       R     G     B\n");
  270.     printf("-----------    ------ ------ ------   ----- ----- -----\n");
  271.     for (t = 1000; t <= 10000; t+= 500) {
  272.  
  273.         bbTemp = t;
  274.         spectrum_to_xyz(bb_spectrum, &x, &y, &z);
  275.         xyz_to_rgb(cs, x, y, z, &r, &g, &b);
  276.         printf("  %5.0f K      %.4f %.4f %.4f   ", t, x, y, z);
  277.         if (constrain_rgb(cs, &x, &y, &z, &r, &g, &b)) {
  278.             printf("%.3f %.3f %.3f (Approximation)\n", r, g, b);
  279.         } else {
  280.             printf("%.3f %.3f %.3f\n", r, g, b);
  281.         }
  282.     }
  283.     return 0;
  284. }
  285.