home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Online Praxis 1996 April
/
OP4_96.ISO
/
surfen
/
astro
/
earthv
/
specrend.c
< prev
next >
Wrap
C/C++ Source or Header
|
1996-08-16
|
12KB
|
285 lines
/*
Colour Rendering of Spectra
by John Walker, kelvin@fourmilab.ch
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;
}