API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Oct 17th, 2019
247
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 15.26 KB | None | 0 0
raw download clone embed print report
  1. /*
  2.                 Colour Rendering of Spectra
  3.  
  4.                        by John Walker
  5.                   http://www.fourmilab.ch/
  6.  
  7.                  Last updated: March 9, 2003
  8.  
  9.            This program is in the public domain.
  10.  
  11.     For complete information about the techniques employed in
  12.     this program, see the World-Wide Web document:
  13.  
  14.              http://www.fourmilab.ch/documents/specrend/
  15.  
  16.     The xyz_to_rgb() function, which was wrong in the original
  17.     version of this program, was corrected by:
  18.  
  19.             Andrew J. S. Hamilton 21 May 1999
  20.             Andrew.Hamilton@Colorado.EDU
  21.             http://casa.colorado.edu/~ajsh/
  22.  
  23.     who also added the gamma correction facilities and
  24.     modified constrain_rgb() to work by desaturating the
  25.     colour by adding white.
  26.  
  27.     A program which uses these functions to plot CIE
  28.     "tongue" diagrams called "ppmcie" is included in
  29.     the Netpbm graphics toolkit:
  30.         http://netpbm.sourceforge.net/
  31.     (The program was called cietoppm in earlier
  32.     versions of Netpbm.)
  33.  
  34. */
  35.  
  36. #include <stdio.h>
  37. #include <math.h>
  38.  
  39. /* A colour system is defined by the CIE x and y coordinates of
  40.    its three primary illuminants and the x and y coordinates of
  41.    the white point. */
  42.  
  43. struct colourSystem {
  44.     char *name;                     /* Colour system name */
  45.     double xRed, yRed,              /* Red x, y */
  46.            xGreen, yGreen,          /* Green x, y */
  47.            xBlue, yBlue,            /* Blue x, y */
  48.            xWhite, yWhite,          /* White point x, y */
  49.            gamma;                   /* Gamma correction for system */
  50. };
  51.  
  52. /* White point chromaticities. */
  53.  
  54. #define IlluminantC     0.3101, 0.3162          /* For NTSC television */
  55. #define IlluminantD65   0.3127, 0.3291          /* For EBU and SMPTE */
  56. #define IlluminantE     0.33333333, 0.33333333  /* CIE equal-energy illuminant */
  57.  
  58. /*  Gamma of nonlinear correction.
  59.  
  60.     See Charles Poynton's ColorFAQ Item 45 and GammaFAQ Item 6 at:
  61.  
  62.        http://www.poynton.com/ColorFAQ.html
  63.        http://www.poynton.com/GammaFAQ.html
  64.  
  65. */
  66.  
  67. #define GAMMA_REC709    0               /* Rec. 709 */
  68.  
  69. static struct colourSystem
  70.                   /* Name                  xRed    yRed    xGreen  yGreen  xBlue  yBlue    White point        Gamma   */
  71.     NTSCsystem  =  { "NTSC",               0.67,   0.33,   0.21,   0.71,   0.14,   0.08,   IlluminantC,    GAMMA_REC709 },
  72.     EBUsystem   =  { "EBU (PAL/SECAM)",    0.64,   0.33,   0.29,   0.60,   0.15,   0.06,   IlluminantD65,  GAMMA_REC709 },
  73.     SMPTEsystem =  { "SMPTE",              0.630,  0.340,  0.310,  0.595,  0.155,  0.070,  IlluminantD65,  GAMMA_REC709 },
  74.     HDTVsystem  =  { "HDTV",               0.670,  0.330,  0.210,  0.710,  0.150,  0.060,  IlluminantD65,  GAMMA_REC709 },
  75.     CIEsystem   =  { "CIE",                0.7355, 0.2645, 0.2658, 0.7243, 0.1669, 0.0085, IlluminantE,    GAMMA_REC709 },
  76.     Rec709system = { "CIE REC 709",        0.64,   0.33,   0.30,   0.60,   0.15,   0.06,   IlluminantD65,  GAMMA_REC709 };
  77.  
  78. /*                          UPVP_TO_XY
  79.  
  80.     Given 1976 coordinates u', v', determine 1931 chromaticities x, y
  81.  
  82. */
  83.  
  84. void upvp_to_xy(double up, double vp, double *xc, double *yc)
  85. {
  86.     *xc = (9 * up) / ((6 * up) - (16 * vp) + 12);
  87.     *yc = (4 * vp) / ((6 * up) - (16 * vp) + 12);
  88. }
  89.  
  90. /*                          XY_TO_UPVP
  91.  
  92.     Given 1931 chromaticities x, y, determine 1976 coordinates u', v'
  93.  
  94. */
  95.  
  96. void xy_to_upvp(double xc, double yc, double *up, double *vp)
  97. {
  98.     *up = (4 * xc) / ((-2 * xc) + (12 * yc) + 3);
  99.     *vp = (9 * yc) / ((-2 * xc) + (12 * yc) + 3);
  100. }
  101.  
  102. /*                             XYZ_TO_RGB
  103.  
  104.     Given an additive tricolour system CS, defined by the CIE x
  105.     and y chromaticities of its three primaries (z is derived
  106.     trivially as 1-(x+y)), and a desired chromaticity (XC, YC,
  107.     ZC) in CIE space, determine the contribution of each
  108.     primary in a linear combination which sums to the desired
  109.     chromaticity.  If the  requested chromaticity falls outside
  110.     the Maxwell  triangle (colour gamut) formed by the three
  111.     primaries, one of the r, g, or b weights will be negative.
  112.  
  113.     Caller can use constrain_rgb() to desaturate an
  114.     outside-gamut colour to the closest representation within
  115.     the available gamut and/or norm_rgb to normalise the RGB
  116.     components so the largest nonzero component has value 1.
  117.  
  118. */
  119.  
  120. void xyz_to_rgb(struct colourSystem *cs,
  121.                 double xc, double yc, double zc,
  122.                 double *r, double *g, double *b)
  123. {
  124.     double xr, yr, zr, xg, yg, zg, xb, yb, zb;
  125.     double xw, yw, zw;
  126.     double rx, ry, rz, gx, gy, gz, bx, by, bz;
  127.     double rw, gw, bw;
  128.  
  129.     xr = cs->xRed;    yr = cs->yRed;    zr = 1 - (xr + yr);
  130.     xg = cs->xGreen;  yg = cs->yGreen;  zg = 1 - (xg + yg);
  131.     xb = cs->xBlue;   yb = cs->yBlue;   zb = 1 - (xb + yb);
  132.  
  133.     xw = cs->xWhite;  yw = cs->yWhite;  zw = 1 - (xw + yw);
  134.  
  135.     /* xyz -> rgb matrix, before scaling to white. */
  136.  
  137.     rx = (yg * zb) - (yb * zg);  ry = (xb * zg) - (xg * zb);  rz = (xg * yb) - (xb * yg);
  138.     gx = (yb * zr) - (yr * zb);  gy = (xr * zb) - (xb * zr);  gz = (xb * yr) - (xr * yb);
  139.     bx = (yr * zg) - (yg * zr);  by = (xg * zr) - (xr * zg);  bz = (xr * yg) - (xg * yr);
  140.  
  141.     /* White scaling factors.
  142.        Dividing by yw scales the white luminance to unity, as conventional. */
  143.  
  144.     rw = ((rx * xw) + (ry * yw) + (rz * zw)) / yw;
  145.     gw = ((gx * xw) + (gy * yw) + (gz * zw)) / yw;
  146.     bw = ((bx * xw) + (by * yw) + (bz * zw)) / yw;
  147.  
  148.     /* xyz -> rgb matrix, correctly scaled to white. */
  149.  
  150.     rx = rx / rw;  ry = ry / rw;  rz = rz / rw;
  151.     gx = gx / gw;  gy = gy / gw;  gz = gz / gw;
  152.     bx = bx / bw;  by = by / bw;  bz = bz / bw;
  153.  
  154.     /* rgb of the desired point */
  155.  
  156.     *r = (rx * xc) + (ry * yc) + (rz * zc);
  157.     *g = (gx * xc) + (gy * yc) + (gz * zc);
  158.     *b = (bx * xc) + (by * yc) + (bz * zc);
  159. }
  160.  
  161. /*                            INSIDE_GAMUT
  162.  
  163.      Test whether a requested colour is within the gamut
  164.      achievable with the primaries of the current colour
  165.      system.  This amounts simply to testing whether all the
  166.      primary weights are non-negative. */
  167.  
  168. int inside_gamut(double r, double g, double b)
  169. {
  170.     return (r >= 0) && (g >= 0) && (b >= 0);
  171. }
  172.  
  173. /*                          CONSTRAIN_RGB
  174.  
  175.     If the requested RGB shade contains a negative weight for
  176.     one of the primaries, it lies outside the colour gamut
  177.     accessible from the given triple of primaries.  Desaturate
  178.     it by adding white, equal quantities of R, G, and B, enough
  179.     to make RGB all positive.  The function returns 1 if the
  180.     components were modified, zero otherwise.
  181.  
  182. */
  183.  
  184. int constrain_rgb(double *r, double *g, double *b)
  185. {
  186.     double w;
  187.  
  188.     /* Amount of white needed is w = - min(0, *r, *g, *b) */
  189.  
  190.     w = (0 < *r) ? 0 : *r;
  191.     w = (w < *g) ? w : *g;
  192.     w = (w < *b) ? w : *b;
  193.     w = -w;
  194.  
  195.     /* Add just enough white to make r, g, b all positive. */
  196.  
  197.     if (w > 0) {
  198.         *r += w;  *g += w; *b += w;
  199.         return 1;                     /* Colour modified to fit RGB gamut */
  200.     }
  201.  
  202.     return 0;                         /* Colour within RGB gamut */
  203. }
  204.  
  205. /*                          GAMMA_CORRECT_RGB
  206.  
  207.     Transform linear RGB values to nonlinear RGB values. Rec.
  208.     709 is ITU-R Recommendation BT. 709 (1990) ``Basic
  209.     Parameter Values for the HDTV Standard for the Studio and
  210.     for International Programme Exchange'', formerly CCIR Rec.
  211.     709. For details see
  212.  
  213.        http://www.poynton.com/ColorFAQ.html
  214.        http://www.poynton.com/GammaFAQ.html
  215. */
  216.  
  217. void gamma_correct(const struct colourSystem *cs, double *c)
  218. {
  219.     double gamma;
  220.  
  221.     gamma = cs->gamma;
  222.  
  223.     if (gamma == GAMMA_REC709) {
  224.         /* Rec. 709 gamma correction. */
  225.         double cc = 0.018;
  226.  
  227.         if (*c < cc) {
  228.             *c *= ((1.099 * pow(cc, 0.45)) - 0.099) / cc;
  229.         } else {
  230.             *c = (1.099 * pow(*c, 0.45)) - 0.099;
  231.         }
  232.     } else {
  233.         /* Nonlinear colour = (Linear colour)^(1/gamma) */
  234.         *c = pow(*c, 1.0 / gamma);
  235.     }
  236. }
  237.  
  238. void gamma_correct_rgb(const struct colourSystem *cs, double *r, double *g, double *b)
  239. {
  240.     gamma_correct(cs, r);
  241.     gamma_correct(cs, g);
  242.     gamma_correct(cs, b);
  243. }
  244.  
  245. /*                          NORM_RGB
  246.  
  247.     Normalise RGB components so the most intense (unless all
  248.     are zero) has a value of 1.
  249.  
  250. */
  251.  
  252. void norm_rgb(double *r, double *g, double *b)
  253. {
  254. #define Max(a, b)   (((a) > (b)) ? (a) : (b))
  255.     double greatest = Max(*r, Max(*g, *b));
  256.  
  257.     if (greatest > 0) {
  258.         *r /= greatest;
  259.         *g /= greatest;
  260.         *b /= greatest;
  261.     }
  262. #undef Max
  263. }
  264.  
  265. /*                          SPECTRUM_TO_XYZ
  266.  
  267.     Calculate the CIE X, Y, and Z coordinates corresponding to
  268.     a light source with spectral distribution given by  the
  269.     function SPEC_INTENS, which is called with a series of
  270.     wavelengths between 380 and 780 nm (the argument is
  271.     expressed in meters), which returns emittance at  that
  272.     wavelength in arbitrary units.  The chromaticity
  273.     coordinates of the spectrum are returned in the x, y, and z
  274.     arguments which respect the identity:
  275.  
  276.             x + y + z = 1.
  277. */
  278.  
  279. void spectrum_to_xyz(double (*spec_intens)(double wavelength),
  280.                      double *x, double *y, double *z)
  281. {
  282.     int i;
  283.     double lambda, X = 0, Y = 0, Z = 0, XYZ;
  284.  
  285.     /* CIE colour matching functions xBar, yBar, and zBar for
  286.        wavelengths from 380 through 780 nanometers, every 5
  287.        nanometers.  For a wavelength lambda in this range:
  288.  
  289.             cie_colour_match[(lambda - 380) / 5][0] = xBar
  290.             cie_colour_match[(lambda - 380) / 5][1] = yBar
  291.             cie_colour_match[(lambda - 380) / 5][2] = zBar
  292.  
  293.         To save memory, this table can be declared as floats
  294.         rather than doubles; (IEEE) float has enough
  295.         significant bits to represent the values. It's declared
  296.         as a double here to avoid warnings about "conversion
  297.         between floating-point types" from certain persnickety
  298.         compilers. */
  299.  
  300.     static double cie_colour_match[81][3] = {
  301.         {0.0014,0.0000,0.0065}, {0.0022,0.0001,0.0105}, {0.0042,0.0001,0.0201},
  302.         {0.0076,0.0002,0.0362}, {0.0143,0.0004,0.0679}, {0.0232,0.0006,0.1102},
  303.         {0.0435,0.0012,0.2074}, {0.0776,0.0022,0.3713}, {0.1344,0.0040,0.6456},
  304.         {0.2148,0.0073,1.0391}, {0.2839,0.0116,1.3856}, {0.3285,0.0168,1.6230},
  305.         {0.3483,0.0230,1.7471}, {0.3481,0.0298,1.7826}, {0.3362,0.0380,1.7721},
  306.         {0.3187,0.0480,1.7441}, {0.2908,0.0600,1.6692}, {0.2511,0.0739,1.5281},
  307.         {0.1954,0.0910,1.2876}, {0.1421,0.1126,1.0419}, {0.0956,0.1390,0.8130},
  308.         {0.0580,0.1693,0.6162}, {0.0320,0.2080,0.4652}, {0.0147,0.2586,0.3533},
  309.         {0.0049,0.3230,0.2720}, {0.0024,0.4073,0.2123}, {0.0093,0.5030,0.1582},
  310.         {0.0291,0.6082,0.1117}, {0.0633,0.7100,0.0782}, {0.1096,0.7932,0.0573},
  311.         {0.1655,0.8620,0.0422}, {0.2257,0.9149,0.0298}, {0.2904,0.9540,0.0203},
  312.         {0.3597,0.9803,0.0134}, {0.4334,0.9950,0.0087}, {0.5121,1.0000,0.0057},
  313.         {0.5945,0.9950,0.0039}, {0.6784,0.9786,0.0027}, {0.7621,0.9520,0.0021},
  314.         {0.8425,0.9154,0.0018}, {0.9163,0.8700,0.0017}, {0.9786,0.8163,0.0014},
  315.         {1.0263,0.7570,0.0011}, {1.0567,0.6949,0.0010}, {1.0622,0.6310,0.0008},
  316.         {1.0456,0.5668,0.0006}, {1.0026,0.5030,0.0003}, {0.9384,0.4412,0.0002},
  317.         {0.8544,0.3810,0.0002}, {0.7514,0.3210,0.0001}, {0.6424,0.2650,0.0000},
  318.         {0.5419,0.2170,0.0000}, {0.4479,0.1750,0.0000}, {0.3608,0.1382,0.0000},
  319.         {0.2835,0.1070,0.0000}, {0.2187,0.0816,0.0000}, {0.1649,0.0610,0.0000},
  320.         {0.1212,0.0446,0.0000}, {0.0874,0.0320,0.0000}, {0.0636,0.0232,0.0000},
  321.         {0.0468,0.0170,0.0000}, {0.0329,0.0119,0.0000}, {0.0227,0.0082,0.0000},
  322.         {0.0158,0.0057,0.0000}, {0.0114,0.0041,0.0000}, {0.0081,0.0029,0.0000},
  323.         {0.0058,0.0021,0.0000}, {0.0041,0.0015,0.0000}, {0.0029,0.0010,0.0000},
  324.         {0.0020,0.0007,0.0000}, {0.0014,0.0005,0.0000}, {0.0010,0.0004,0.0000},
  325.         {0.0007,0.0002,0.0000}, {0.0005,0.0002,0.0000}, {0.0003,0.0001,0.0000},
  326.         {0.0002,0.0001,0.0000}, {0.0002,0.0001,0.0000}, {0.0001,0.0000,0.0000},
  327.         {0.0001,0.0000,0.0000}, {0.0001,0.0000,0.0000}, {0.0000,0.0000,0.0000}
  328.     };
  329.  
  330.     for (i = 0, lambda = 380; lambda < 780.1; i++, lambda += 5) {
  331.         double Me;
  332.  
  333.         Me = (*spec_intens)(lambda);
  334.         X += Me * cie_colour_match[i][0];
  335.         Y += Me * cie_colour_match[i][1];
  336.         Z += Me * cie_colour_match[i][2];
  337.     }
  338.     XYZ = (X + Y + Z);
  339.     *x = X / XYZ;
  340.     *y = Y / XYZ;
  341.     *z = Z / XYZ;
  342. }
  343.  
  344. /*                            BB_SPECTRUM
  345.  
  346.     Calculate, by Planck's radiation law, the emittance of a black body
  347.     of temperature bbTemp at the given wavelength (in metres).  */
  348.  
  349. double bbTemp = 5000;                 /* Hidden temperature argument
  350.                                          to BB_SPECTRUM. */
  351. double bb_spectrum(double wavelength)
  352. {
  353.     double wlm = wavelength * 1e-9;   /* Wavelength in meters */
  354.  
  355.     return (3.74183e-16 * pow(wlm, -5.0)) /
  356.            (exp(1.4388e-2 / (wlm * bbTemp)) - 1.0);
  357. }
  358.  
  359. /*  Built-in test program which displays the x, y, and Z and RGB
  360.     values for black body spectra from 1000 to 10000 degrees kelvin.
  361.     When run, this program should produce the following output:
  362.  
  363.     Temperature       x      y      z       R     G     B
  364.     -----------    ------ ------ ------   ----- ----- -----
  365.        1000 K      0.6528 0.3444 0.0028   1.000 0.007 0.000 (Approximation)
  366.        1500 K      0.5857 0.3931 0.0212   1.000 0.126 0.000 (Approximation)
  367.        2000 K      0.5267 0.4133 0.0600   1.000 0.234 0.010
  368.        2500 K      0.4770 0.4137 0.1093   1.000 0.349 0.067
  369.        3000 K      0.4369 0.4041 0.1590   1.000 0.454 0.151
  370.        3500 K      0.4053 0.3907 0.2040   1.000 0.549 0.254
  371.        4000 K      0.3805 0.3768 0.2428   1.000 0.635 0.370
  372.        4500 K      0.3608 0.3636 0.2756   1.000 0.710 0.493
  373.        5000 K      0.3451 0.3516 0.3032   1.000 0.778 0.620
  374.        5500 K      0.3325 0.3411 0.3265   1.000 0.837 0.746
  375.        6000 K      0.3221 0.3318 0.3461   1.000 0.890 0.869
  376.        6500 K      0.3135 0.3237 0.3628   1.000 0.937 0.988
  377.        7000 K      0.3064 0.3166 0.3770   0.907 0.888 1.000
  378.        7500 K      0.3004 0.3103 0.3893   0.827 0.839 1.000
  379.        8000 K      0.2952 0.3048 0.4000   0.762 0.800 1.000
  380.        8500 K      0.2908 0.3000 0.4093   0.711 0.766 1.000
  381.        9000 K      0.2869 0.2956 0.4174   0.668 0.738 1.000
  382.        9500 K      0.2836 0.2918 0.4246   0.632 0.714 1.000
  383.       10000 K      0.2807 0.2884 0.4310   0.602 0.693 1.000
  384. */
  385.  
  386. int main()
  387. {
  388.     double t, x, y, z, r, g, b;
  389.     struct colourSystem *cs = &SMPTEsystem;
  390.  
  391.     printf("Temperature       x      y      z       R     G     B\n");
  392.     printf("-----------    ------ ------ ------   ----- ----- -----\n");
  393.  
  394.     for (t = 1000; t <= 10000; t+= 500) {
  395.         bbTemp = t;
  396.         spectrum_to_xyz(bb_spectrum, &x, &y, &z);
  397.         xyz_to_rgb(cs, x, y, z, &r, &g, &b);
  398.         printf("  %5.0f K      %.4f %.4f %.4f   ", t, x, y, z);
  399.         if (constrain_rgb(&r, &g, &b)) {
  400.             norm_rgb(&r, &g, &b);
  401.             printf("%.3f %.3f %.3f (Approximation)\n", r, g, b);
  402.         } else {
  403.             norm_rgb(&r, &g, &b);
  404.             printf("%.3f %.3f %.3f\n", r, g, b);
  405.         }
  406.     }
  407.     return 0;
  408. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!