Re: темпеpатypа->цвет

From
Alex Ivchenko ()
To
Peter Pike's
Date
2002-07-04T00:07:30Z
Area
RU.ALGORITHMS
From: Alex Ivchenko <alexx@ueidaq.com>

Peter,

> Hyжно наглядно пpеобpазовать темпеpатypнyю шкалy в цветовyю (256 цветов).
> Сабжевый алгоpитм обсyждался в эхе, но y меня все отпypжилось.
> Не подскажете, как это можно сделать?

You can derive calculations from this one:

--A

/*
                     Colour Rendering of Spectra

                            by John Walker
                       http://www.fourmilab.ch/

                This program is in the public domain.

    For complete information about the techniques employed in this
    program, see the World-Wide Web document:

             http://www.fourmilab.ch/documents/specrend/

*/

#include <stdio.h>
#include <math.h>

#define CLAMP(v,l,h)    ((v)<(l) ? (l) : (v) > (h) ? (h) : (v))

/* A colour system is defined by the CIE x and y coordinates of its
   three primary illuminants and the x and y coordinates of the white
   point. */

struct colourSystem {
    double xRed, yRed,
           xGreen, yGreen,
           xBlue, yBlue,
           xWhite, yWhite;
};

/* White point chromaticities. */

#define IlluminantC     0.3101, 0.3162  /* For NTSC television */
#define IlluminantD65   0.3127, 0.3291  /* For EBU and SMPTE */

static struct colourSystem
                  /* xRed   yRed  xGreen yGreen  xBlue  yBlue  White point */
    NTSCsystem  =  { 0.67,  0.33,  0.21,  0.71,  0.14,  0.08,  IlluminantC   },
    EBUsystem   =  { 0.64,  0.33,  0.29,  0.60,  0.15,  0.06,  IlluminantD65 },
    SMPTEsystem =  { 0.630, 0.340, 0.310, 0.595, 0.155, 0.070, IlluminantD65 };

/*                             XYZ_TO_RGB

    Given an additive tricolour system CS, defined by the  CIE  x  and  y
    chromaticities  of  its  three  primaries (z is derived trivially as
    1-(x+y)), and a desired chromaticity (XC,  YC,  ZC)  in  CIE  space,
    determine  the  contribution of each primary in a linear combination
    which  sums  to  the  desired  chromaticity.    If   the   requested
    chromaticity falls outside the Maxwell triangle (colour gamut) formed
    by the three primaries, one of the  r,  g,  or  b  weights  will  be
    negative.   Use  inside_gamut()  to  test  for  a  valid  colour  and
    constrain_rgb() to desaturate an outside-gamut colour to the  closest
    representation within the available gamut. */


void xyz_to_rgb(struct colourSystem *cs,
                double xc, double yc, double zc,
                double *r, double *g, double *b)
{
    double xr, yr, zr, xg, yg, zg, xb, yb, zb, d;

    xr = cs->xRed;    yr = cs->yRed;    zr = 1 - (xr + yr);
    xg = cs->xGreen;  yg = cs->yGreen;  zg = 1 - (xg + yg);
    xb = cs->xBlue;   yb = cs->yBlue;   zb = 1 - (xb + yb);
    d = xr*yg*zb - xg*yr*zb - xr*yb*zg + xb*yr*zg + xg*yb*zr - xb*yg*zr;

    *r = (-xg*yc*zb + xc*yg*zb + xg*yb*zc - xb*yg*zc - xc*yb*zg + xb*yc*zg) / d;

    *g = (xr*yc*zb - xc*yr*zb - xr*yb*zc + xb*yr*zc + xc*yb*zr - xb*yc*zr) / d;

    *b = (xr*yg*zc - xg*yr*zc - xr*yc*zg + xc*yr*zg + xg*yc*zr - xc*yg*zr) / d;
}

/*                            INSIDE_GAMUT

    Test  whether  a requested colour is within the gamut achievable with
    the primaries of the current colour system.  This amounts  simply  to
    testing whether all the primary weights are non-negative. */

int inside_gamut(double r, double g, double b)
{
    return r >= 0 && g >= 0 && b >= 0;
}

/*                           CONSTRAIN_RGB

    If  the  requested RGB shade contains a negative weight for one of
    the primaries, it lies outside the colour gamut accessible from the
    given triple of primaries.  Desaturate it by mixing with the white
    point of the colour system so as to reduce  the  primary  with  the
    negative weight to  zero.   This  is  equivalent  to  finding  the
    intersection  on the CIE diagram of a line drawn between the white
    point and the  requested  colour  with  the  edge  of  the  Maxwell
    triangle formed by the three primaries. */

int constrain_rgb(struct colourSystem *cs,
                  double *x, double *y, double *z,
                  double *r, double *g, double *b)
{
    /* Is the contribution of one of the primaries negative ? */

    if (!inside_gamut(*r, *g, *b)) {
        double par, wr, wg, wb;

        /* Yes.  Find the RGB mixing weights of the white point (we
           assume the white point is in the gamut!). */

        xyz_to_rgb(cs,
                   cs->xWhite, cs->yWhite, 1 - (cs->xWhite + cs->yWhite),
                   &wr, &wg, &wb);

        /* Find the primary with negative weight and calculate the
           parameter of the point on the vector from the white point
           to the original requested colour in RGB space. */

        if (*r < *g && *r < *b) {
            par = wr / (wr - *r);
        } else if (*g < *r && *g < *b) {
            par = wg / (wg - *g);
        } else {
            par = wb / (wb - *b);
        }

        /* Since XYZ space is a linear transformation of RGB space, we
           can find the XYZ space coordinates of the point where the
           edge of the gamut intersects the vector from the white point
           to the original colour by multiplying the parameter in RGB
           space by the difference vector in XYZ space. */

        *x = CLAMP(cs->xWhite + par * (*x - cs->xWhite), 0, 1);
        *y = CLAMP(cs->yWhite + par * (*y - cs->yWhite), 0, 1);
        *z = CLAMP(1 - (*x + *y), 0, 1);

        /* Now finally calculate the gamut-constrained RGB weights. */

        *r = CLAMP(wr + par * (*r - wr), 0, 1);
        *g = CLAMP(wg + par * (*g - wg), 0, 1);
        *b = CLAMP(wb + par * (*b - wb), 0, 1);
        return 1;                     /* Colour modified to fit RGB gamut */
    }
    return 0;                         /* Colour within RGB gamut */
}

/*                          SPECTRUM_TO_XYZ

    Calculate the CIE X, Y, and Z coordinates corresponding to a light
    source  with  spectral  distribution   given   by   the   function
    SPEC_INTENS,  which is called with a series of wavelengths between
    380 and 780 nm  (the  argument  is  expressed  in  meters),  which
    returns  emittance  at  that  wavelength  in arbitrary units.  The
    chromaticity coordinates of the spectrum are returned in the x, y,
    and z arguments which respect the identity:

            x + y + z = 1.
*/

void spectrum_to_xyz(double (*spec_intens)(double wavelength),
                     double *x, double *y, double *z)
{
    int i;
    double lambda, X = 0, Y = 0, Z = 0, XYZ;

    /* CIE colour matching functions xBar, yBar, and zBar for
       wavelengths from 380 through 780 nanometers, every 5
       nanometers.  For a wavelength lambda in this range:

            cie_colour_match[(lambda - 380) / 5][0] = xBar
            cie_colour_match[(lambda - 380) / 5][1] = yBar
            cie_colour_match[(lambda - 380) / 5][2] = zBar

       To  save  memory,  this  table can be declared as floats rather
       than doubles; (IEEE)  float  has  enough  significant  bits  to
       represent  the values.  It's declared as a double here to avoid
       warnings about "conversion between floating-point  types"  from
       certain persnickety compilers. */

    static double cie_colour_match[81][3] = {
        {0.0014,0.0000,0.0065}, {0.0022,0.0001,0.0105}, {0.0042,0.0001,0.0201},
        {0.0076,0.0002,0.0362}, {0.0143,0.0004,0.0679}, {0.0232,0.0006,0.1102},
        {0.0435,0.0012,0.2074}, {0.0776,0.0022,0.3713}, {0.1344,0.0040,0.6456},
        {0.2148,0.0073,1.0391}, {0.2839,0.0116,1.3856}, {0.3285,0.0168,1.6230},
        {0.3483,0.0230,1.7471}, {0.3481,0.0298,1.7826}, {0.3362,0.0380,1.7721},
        {0.3187,0.0480,1.7441}, {0.2908,0.0600,1.6692}, {0.2511,0.0739,1.5281},
        {0.1954,0.0910,1.2876}, {0.1421,0.1126,1.0419}, {0.0956,0.1390,0.8130},
        {0.0580,0.1693,0.6162}, {0.0320,0.2080,0.4652}, {0.0147,0.2586,0.3533},
        {0.0049,0.3230,0.2720}, {0.0024,0.4073,0.2123}, {0.0093,0.5030,0.1582},
        {0.0291,0.6082,0.1117}, {0.0633,0.7100,0.0782}, {0.1096,0.7932,0.0573},
        {0.1655,0.8620,0.0422}, {0.2257,0.9149,0.0298}, {0.2904,0.9540,0.0203},
        {0.3597,0.9803,0.0134}, {0.4334,0.9950,0.0087}, {0.5121,1.0000,0.0057},
        {0.5945,0.9950,0.0039}, {0.6784,0.9786,0.0027}, {0.7621,0.9520,0.0021},
        {0.8425,0.9154,0.0018}, {0.9163,0.8700,0.0017}, {0.9786,0.8163,0.0014},
        {1.0263,0.7570,0.0011}, {1.0567,0.6949,0.0010}, {1.0622,0.6310,0.0008},
        {1.0456,0.5668,0.0006}, {1.0026,0.5030,0.0003}, {0.9384,0.4412,0.0002},
        {0.8544,0.3810,0.0002}, {0.7514,0.3210,0.0001}, {0.6424,0.2650,0.0000},
        {0.5419,0.2170,0.0000}, {0.4479,0.1750,0.0000}, {0.3608,0.1382,0.0000},
        {0.2835,0.1070,0.0000}, {0.2187,0.0816,0.0000}, {0.1649,0.0610,0.0000},
        {0.1212,0.0446,0.0000}, {0.0874,0.0320,0.0000}, {0.0636,0.0232,0.0000},
        {0.0468,0.0170,0.0000}, {0.0329,0.0119,0.0000}, {0.0227,0.0082,0.0000},
        {0.0158,0.0057,0.0000}, {0.0114,0.0041,0.0000}, {0.0081,0.0029,0.0000},
        {0.0058,0.0021,0.0000}, {0.0041,0.0015,0.0000}, {0.0029,0.0010,0.0000},
        {0.0020,0.0007,0.0000}, {0.0014,0.0005,0.0000}, {0.0010,0.0004,0.0000},
        {0.0007,0.0002,0.0000}, {0.0005,0.0002,0.0000}, {0.0003,0.0001,0.0000},
        {0.0002,0.0001,0.0000}, {0.0002,0.0001,0.0000}, {0.0001,0.0000,0.0000},
        {0.0001,0.0000,0.0000}, {0.0001,0.0000,0.0000}, {0.0000,0.0000,0.0000}
    };

    for (i = 0, lambda = 380; lambda < 780.1; i++, lambda += 5) {
        double Me;

        Me = (*spec_intens)(lambda);
        X += Me * cie_colour_match[i][0];
        Y += Me * cie_colour_match[i][1];
        Z += Me * cie_colour_match[i][2];
    }
    XYZ = (X + Y + Z);
    *x = X / XYZ;
    *y = Y / XYZ;
    *z = Z / XYZ;
}

/*                            BB_SPECTRUM

    Calculate, by Planck's radiation law, the emittance of a black body
    of temperature bbTemp at the given wavelength (in meters).  */

double bbTemp = 5000;                 /* Hidden temperature argument
                                         to BB_SPECTRUM. */
double bb_spectrum(double wavelength)
{
    double wlm = wavelength * 1e-9;   /* Wavelength in meters */

    return (3.74183e-16 * pow(wlm, -5.0)) /
           (exp(1.4388e-2 / (wlm * bbTemp)) - 1.0);
}

/*  Built-in test program which displays the x, y, and Z and RGB
    values for black body spectra from 1000 to 10000 degrees kelvin.
    When run, this program should produce the following output:

    Temperature       x      y      z       R     G     B
    -----------    ------ ------ ------   ----- ----- -----
       1000 K      0.6528 0.3444 0.0028   0.988 0.012 0.000 (Approximation)
       1500 K      0.5857 0.3931 0.0212   0.808 0.192 0.000 (Approximation)
       2000 K      0.5267 0.4133 0.0600   0.684 0.302 0.014
       2500 K      0.4770 0.4137 0.1093   0.558 0.368 0.074
       3000 K      0.4369 0.4041 0.1590   0.464 0.398 0.138
       3500 K      0.4053 0.3907 0.2040   0.394 0.408 0.198
       4000 K      0.3805 0.3768 0.2428   0.341 0.409 0.250
       4500 K      0.3608 0.3636 0.2756   0.301 0.404 0.294
       5000 K      0.3451 0.3516 0.3032   0.271 0.397 0.332
       5500 K      0.3325 0.3411 0.3265   0.246 0.390 0.364
       6000 K      0.3221 0.3318 0.3461   0.227 0.382 0.391
       6500 K      0.3135 0.3237 0.3628   0.212 0.374 0.414
       7000 K      0.3064 0.3166 0.3770   0.199 0.367 0.434
       7500 K      0.3004 0.3103 0.3893   0.188 0.361 0.451
       8000 K      0.2952 0.3048 0.4000   0.179 0.355 0.466
       8500 K      0.2908 0.3000 0.4093   0.172 0.350 0.479
       9000 K      0.2869 0.2956 0.4174   0.165 0.345 0.490
       9500 K      0.2836 0.2918 0.4246   0.160 0.340 0.500
      10000 K      0.2807 0.2884 0.4310   0.155 0.336 0.509
*/

int main()
{ 
    double t, x, y, z, r, g, b;
    struct colourSystem *cs = &SMPTEsystem;

    printf("Temperature       x      y      z       R     G     B\n");
    printf("-----------    ------ ------ ------   ----- ----- -----\n");
    for (t = 1000; t <= 10000; t+= 500) {

        bbTemp = t;
        spectrum_to_xyz(bb_spectrum, &x, &y, &z);
        xyz_to_rgb(cs, x, y, z, &r, &g, &b);
        printf("  %5.0f K      %.4f %.4f %.4f   ", t, x, y, z);
        if (constrain_rgb(cs, &x, &y, &z, &r, &g, &b)) {
            printf("%.3f %.3f %.3f (Approximation)\n", r, g, b);
        } else {
            printf("%.3f %.3f %.3f\n", r, g, b);
        }
    }
    return 0;
}
--- ifmail v.2.15dev5
 * Origin: UEI, Inc. (2:5020/400)