Spectrum of Blackbodies

This lesson will provide a very quick introduction to the computation of blackbody's spectrum. It uses a few techniques related to spectrum rendering, color conversion, color space etc. which, due to the complexity of these topics, are studied in other lessons. The reader is assumed to have a good understanding of these techniques which are not explained in this lesson.

Blackbodies refer to objects which have the property to absorb almost all incident radiation and re-emit this energy in the form of light. This is a very simple (not scientific) definition to describe in short what we know as fire, incandescent metals, or even light bulbs. All these objects have the property of emitting light which is the result of a radiative process (as we are not physicists we won't be exploring this complex topic but if you are an expert on thermodynamics we welcome an explanation). One of the most important characteristics of these objects is that the spectrum emitted (or radiated) depends on the object's temperature (you should be familiar with the concept of spectrum which is in short, a wavelength representation of color). When the temperature of a blackbody approaches 1000K (Kelvin degrees) it takes a deep red color. As the temperature increases, it becomes more yellowish (4500K) then almost white (6500K) then becomes pale blue as it reaches 8000 to 9000K. It is important to note that almost all liquid and solid objects start to glow when they reach a temperature of about 800K. Objects emit light below this temperature however it is too weak to be perceivable. 

The Mystery of Light Color

Figure 1: a photograph of various lamps illustrates the effect of colour temperature differences. From left to right: fluorescent light (6500K), incandescent, fluorescent (2644K), fluorescent (3000K).

We are all familiar with light bulbs in which a filament is heated to a temperature at which a fraction of the radiation falls in the visible spectrum. The color of the light depends indeed on the filaments temperature which explains why light bulb color varies from orange to blue (assuming the glass bulb is transparent or white). It also explains the term color temperature which we use in photography and which refers to the color effect we get from using a specific type of light when we take a photograph. Candlelight's temperature varies from 1000 to 2000K. Typical incandescent lamp have a temperature around 2500K. The sun's surface temperature is about 6500K (the sun light color is actually white). Light emitted by objects which temperature is around 2500-3500K is said to be warm (red/orange) while light emitted by objects which temperature is around 6500 is said to be cold (blue/white). It is interesting to know that the peak emittance wavelength for incandescent lamps (around 1160 nm) is actually in the infrared, a part of the spectrum that is not visible to the human eye.

Plank's Law

The power spectral density (in other words how bright is the object) for a given wavelength and temperature can be computed using Plank's law:

$$\begin{array}{ll} I(\lambda,T)=\frac{2hc^2}{\lambda^5}\frac{1}{e^{\frac{hc}{\lambda kT}}-1} \end{array}$$

where h is the Plank constant (xx), c is the speed of light (xx ms-1), k is the Boltzmann constant (xx), λ is the wavelength and T is the temperature of the blackbody. To compute the color of a blackbody for a given temperature, all we need to do is to evaluate Plank's equation from 380 nm to 780 nm (which covers the visible spectrum). We will then convert the resulting blackbody's spectrum to XYZ tristimulus using the color matching function and finally convert the XYZ color to RGB (spectrum and color conversion are described in the lesson on Color in the basic section).

// insert figure of curves different temperatures xx

Our example program will compute this color for a range of temperature varying from 1000K to 10000K. If the XYZ color is not scaled properly, some RGB colors might be greater than 1. Care must be taken to scale the colors down before you save the result to a bitmap image. Another solution is to use a floating point image format (PFM, OpenEXR, Radiance HDR). Some RGB colors can also be negative. It is possible to alter the position of these colors by forcing them to lie in the gamut of the chosen color space however, in our example, for simplicity, we will just clamp negative values.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
// CIE color matching functions static double colorMatchingFunc[][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} }; const double wavelengthMin = 380; const double wavelengthMax = 780; const double wavelengthStep = 5; const unsigned numSteps = (wavelengthMax - wavelengthMin) / wavelengthStep + 1; inline double computePlanck( const double &T, // temperature (Kelvin) const double &lambda) // wavelength (meter) { static const double h = 6.62606896e-34; // Plank constant static const double c = 2.99792458e+8; // Speed of light static const double k = 1.38064880e-23; // Boltzmann constant static const double arg1 = 2 * M_PI * h * c * c; static const double arg2 = ( h * c ) / k; return (arg1 * pow(lambda, -5.0)) / (exp(arg2 / (lambda * T)) - 1.0); } template<typename T> class Color3 { public: Color3() : r(0), g(0), b(0) {} Color3( T rr, T gg, T bb) : r(rr), g(gg), b(bb) {} T r, g, b; }; typedef Color3<double> Color3d; // convert XYZ tristimulus values to RGB Color3d XYZ2RGB(const double &X, const double &Y, const double &Z) { // convert XYZ to sRGB return Color3d( X * 3.2406 + Y * -1.5372 + Z * -0.4986, X * -0.9689 + Y * 1.8758 + Z * 0.0415, X * 0.0557 + Y * -0.2040 + Z * 1.0570 ); } Color3d computeBlackbody(const double &temperature, double &power) { power = 0; double X = 0, Y = 0, Z = 0; for (unsigned k = 0; k < numSteps; k++ ) { // convert to nanometre double lambda = (wavelengthMin + k * wavelengthStep) * 1e-9; double I = computePlanck(temperature, lambda); power += I * wavelengthStep; X += I * colorMatchingFunc[k][0]; Y += I * colorMatchingFunc[k][1]; Z += I * colorMatchingFunc[k][2]; } power /= double( numSteps ); // normalise the result double nor = 1 / std::max( X, std::max( Y, Z ) ); X *= nor, Y *= nor, Z *= nor; return XYZ2RGB(X, Y, Z); } void computeBlackbodyRamp() { unsigned rampWidth = 512, rampHeight = 50; unsigned tempMin = 1000; unsigned tempMax = 10000; Color3d result[ rampWidth ]; double power, temperature; for (unsigned i = 0; i < rampWidth; ++i) { temperature = tempMin + (tempMax - tempMin) * i / double(rampWidth - 1); result[ i ] = computeBlackbody(temperature, power ); } // convert double rgb to unsigned bytes unsigned numBytes = 3 * rampWidth; unsigned char lines[ numBytes ]; unsigned r, g, b; double scale = 1 / 2.6; // arbitrary division to keep values < 1 for (unsigned i = 0, k = 0; i < rampWidth; ++i, k += 3) { r = (unsigned char)(std::max( 0., std::min( 1., result[i].r * scale ) ) * 255 + 0.5); g = (unsigned char)(std::max( 0., std::min( 1., result[i].g * scale ) ) * 255 + 0.5); b = (unsigned char)(std::max( 0., std::min( 1., result[i].b * scale ) ) * 255 + 0.5); lines[ k ] = r, lines[ k + 1 ] = g, lines[ k + 2 ] = b; } // save the result to PPM file std::ofstream ofs("./blackbodyramp.ppm"); ofs << "P6\n" << rampWidth << " " << rampHeight << "\n255\n"; for (unsigned j = 0; j < rampHeight; ++j) { ofs.write((char*)lines, numBytes); } ofs.close(); }

Stefan-Boltzmann's Law

So far we learned how  to convert a blackbody to a color based on its temperature. However in a renderer, it might be useful to measure the amount of energy that is emitted by this blackbody as well. There's a couple of options we can use. Each value we have computed in the spectrum using the Planck radiation formula, represents the amount of energy emitted by the blackbody for a certain wavelength. To find the total amount of energy emitted, all we need to do is sum up all these values. This is a classical numerical integration: the area under the curve/profile of the spectrum we have computed corresponds to the the total power radiated by the blackbody (check the lesson on Numerical Integration in the advanced lessons). Take each one the values from the spectrum, multiplied them by dλ (the wavelength interval which is 5) sum them up and divide this sum by the total number of intervals (81). Note that the result of this integration only accounts for the power radiated within the range of wavelengths corresponding to the visible spectrum (380-780 nm). It is very likely that this number represents a fraction only of the total power emitted over the entire spectrum. What you get is a radiated power (radiant flux or radiant power) expressed in watts (per unit time ? to be confirmed).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Color3d computeBlackbody( const double &temperature, double &power ) { power = 0; ... for (unsigned k = 0; k < numSteps; ++k) { // convert to nanometre double lambda = (wavelengthMin + k * wavelengthStep) * 1e-9; double I = computePlanck(temperature, lambda); power += I * wavelengthStep; ... } power /= double( numSteps ); ... }

The energy radiated by a blackbody can also be computed analytically using the Stefan-Boltzmann's law. Its states that the energy radiated by a blackbody radiator per second per unit area is proportional to the fourth power of the absolute temperature and is given by:

$$\begin{array}{ll} \sigma T^4 \quad j/m^2s\\ \sigma=5.6703*10^-8 \quad watt/m^2K^4 \end{array}$$

Figure 2: in red, the result of the total power emitted by a blackbody function of the temperature using numerical integration of the Plank's formula. In blue, we used the Stefan-Boltzmann's formula to compute the same value. The difference between the two curves is partially due to the fact that we only render the power for a range of wavelengths contained in the visible spectrum (in red) while the Stefan-Boltzmann equation takes into the entire spectrum (in blue). What's important though is that the two curves display an identical shape.

The Planckian locus

When colors emitted by a blackbody for increasing temperature are connected on a chromaticity diagram, they form a curve which is known as the Planckian Locus.

Figure 3: Planckian locus curve plotted on top of the CIE xy 1931 chromaticity diagram

You can use the Planck's law to compute this curve or compute a fast approximation using the following formula (in code):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
Color3d computePlanckLocus( const double &T ) { double xc, yc, arg1 = 1e9 / (T * T * T), arg2 = 1e6 / (T * T), arg3 = 1e3 / T; if ( T >= 1667 && T <= 4000 ) { xc = -0.2661239 * arg1 - 0.2343580 * arg2 + 0.8776956 * arg3 + 0.179910; } else if ( T > 4000 && T <= 25000 ) { xc = -3.0258469 * arg1 + 2.1070379 * arg2 + 0.2226347 * arg3 + 0.240390; } double xc3=xc*xc*xc, xc2=xc*xc; if ( T >= 1667 && T <= 2222 ) { yc = -1.1063814 * xc3 - 1.34811020 * xc2 + 2.18555832 * xc - 0.20219683; } else if ( T > 2222 && T <= 4000 ) { yc = -0.9549476 * xc3 - 1.37418593 * xc2 + 2.09137015 * xc - 0.16748867; } else if ( T > 4000 && T <= 25000 ) { yc = +3.0817580 * xc3 - 5.87338670 * xc2 + 3.75112997 * xc - 0.37001483; } double X = xc / yc; double Y = 1; double Z = (1 - xc - yc) / yc; return XYZ2RGB(X, Y, Z); }

As shown below, the result is not as accurate as the ramp we computed using the spectrum technique (however this can be a matter of interpretation). The result is still quite good and the function is about 66 times faster to compute than the spectrum version. For realtime applications, it might be preferable to use this equation combined with the Stefan-Boltzmann law to evaluate the overall intensity emitted by the blackbody.

Figure 4: a comparaison of the colors computed with the Planckian Locus formula and the spectrum. The spectrum ramp (bottom) seems to give more accurate and realistic colors. However the Planckian locus (top) gives a good approximation and is about 66 faster to compute (on the same machine) than the spectrum version which is an advantage for realtime simulations.

Chapter 1 of 1