## Understanding blackbody radiation for realistic color and light in CG simulations.

**Reading time: 11 mins.**

## Spectrum of Blackbodies

This lesson provides a concise introduction to computing the spectrum of blackbodies. It employs several techniques related to spectrum rendering, color conversion, and color space. Due to the complexity of these topics, they are explored in detail in other lessons. Readers are assumed to possess a solid understanding of these techniques, which are not explained within this lesson.

Blackbodies are objects characterized by their ability to absorb nearly all incident radiation and re-emit this energy as light. This definition, while simplified and not deeply scientific, helps describe phenomena such as fire, incandescent metals, or even light bulbs. These objects emit light as a result of a radiative process. While the intricacies of this process are beyond our scope (as we do not delve into physics here), experts in thermodynamics are welcome to provide further explanation. A key feature of these objects is that the spectrum they emit depends on their temperature. The concept of a spectrum, briefly, refers to a representation of color by wavelength.

As the temperature of a blackbody reaches approximately 1000K (Kelvin), it exhibits a deep red color. With increasing temperature, the color shifts to yellowish at around 4500K, approaches white near 6500K, and transitions to pale blue at temperatures between 8000 and 9000K. It's noteworthy that almost all liquid and solid objects begin to glow when they reach temperatures about 800K. Although objects emit light below this threshold, the emitted light is too weak to be noticeable.

## The Mystery of Light Color

We often encounter light bulbs where a filament is heated to the point that it emits radiation, some of which falls within the visible spectrum. The color of the emitted light directly correlates with the filament's temperature, leading to a variation in light bulb color from orange to blue, assuming the glass bulb is either transparent or white. This phenomenon underpins the concept of **color temperature** used in photography, referring to the color effect obtained from using a specific type of light in photographing subjects. Candlelight, for example, exhibits a temperature ranging from 1000 to 2000K, suggesting a warm glow. Typical incandescent lamps operate around 2500K, emitting a similarly warm light. In contrast, the sun's surface temperature is about 6500K, rendering sunlight as white light.

Light from objects with temperatures around 2500-3500K is considered warm, characterized by red/orange hues, while light from sources around 6500K is deemed cold, with blue/white tones. An interesting point to note is that the peak emittance wavelength for incandescent lamps, around 1160 nm, actually lies in the infrared portion of the spectrum, which is invisible to the human eye. This infrared radiation contributes significantly to the heat emitted by such lamps, beyond the visible light they produce.

## Planck's Law

The power spectral density (essentially, how bright an object is) for a given wavelength and temperature can be computed using **Planck's law**:

where \(h\) is the Planck constant (xx), \(c\) is the speed of light (xx m/s), \(k\) is the Boltzmann constant (xx), \(\lambda\) is the wavelength, and \(T\) is the temperature of the blackbody. To compute the color of a blackbody for a given temperature, we need to evaluate Planck's equation from 380 nm to 780 nm (which covers the visible spectrum). We will then convert the resulting blackbody spectrum to XYZ tristimulus values 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).

Our example program will compute this color for a range of temperatures varying from 1000K to 10000K. If the XYZ color is not scaled properly, some RGB colors might exceed 1. Care must be taken to scale the colors down before saving 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. While it is possible to alter the position of these colors to fit within the gamut of the chosen color space, in our example, for simplicity, we will just clamp negative values.

// 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 ); // Normalize 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 Law

We've learned how to convert a blackbody's temperature into color. In rendering, it's also useful to determine the amount of energy emitted by this blackbody. There are a couple of methods we can employ. The values calculated in the spectrum using the Planck radiation formula each represent the energy emitted by the blackbody at a specific wavelength. To find the total emitted energy, we simply sum these values. This process is a classic example of numerical integration: the area under the curve/profile of the spectrum we've computed corresponds to the total power radiated by the blackbody (refer to The Mathematics of Shading lesson). Multiply each value from the spectrum by \( \Delta\lambda \) (the wavelength interval, which is 5), sum them up, and divide this sum by the total number of intervals (81). Note that this integration result only accounts for the power radiated within the visible spectrum range (380-780 nm) and likely represents only a fraction of the total power emitted across the entire spectrum. The result is the radiated power (also known as **radiant flux** or **radiant power**), expressed in watts (per unit time? to be confirmed).

Color3d computeBlackbody(const double &temperature, double &power) { power = 0; ... for (unsigned k = 0; k < numSteps; ++k) { // convert to nanometers 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 analytically calculated using Stefan-Boltzmann law. It states that the energy radiated per second per unit area (joules per second per square meter, or, watts per square meter) by a blackbody radiator is proportional to the fourth power of the absolute temperature (measured in Kelvin) and is given by:

$$\sigma T^4$$where the constant \(\sigma\), is called the Stefanâ€“Boltzmann constant. It has a value of

$$\sigma = 5.670374419 \times 10^{-8} \, \text{W} \cdot \text{m}^{-2} \cdot \text{K}^{-4}.$$## The Planckian Locus

On a chromaticity diagram, the curve formed by connecting the colors emitted by a blackbody as its temperature increases is known as the Planckian locus.

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

Color3d computePlanckLocus(const double &T) { double x_c, y_c; double arg1 = 1e9 / (T * T * T), arg2 = 1e6 / (T * T), arg3 = 1e3 / T; if (T >= 1667 && T <= 4000) { x_c = -0.2661239 * arg1 - 0.2343580 * arg2 + 0.8776956 * arg3 + 0.179910; } else if (T > 4000 && T <= 25000) { x_c = -3.0258469 * arg1 + 2.1070379 * arg2 + 0.2226347 * arg3 + 0.240390; } double x_c3 = x_c * x_c * x_c, x_c2 = x_c * x_c; if (T >= 1667 && T <= 2222) { y_c = -1.1063814 * x_c3 - 1.34811020 * x_c2 + 2.18555832 * x_c - 0.20219683; } else if (T > 2222 && T <= 4000) { y_c = -0.9549476 * x_c3 - 1.37418593 * x_c2 + 2.09137015 * x_c - 0.16748867; } else if (T > 4000 && T <= 25000) { y_c = +3.0817580 * x_c3 - 5.87338670 * x_c2 + 3.75112997 * x_c - 0.37001483; } double X = x_c / y_c; double Y = 1; double Z = (1 - x_c - y_c) / y_c; return XYZ2RGB(X, Y, Z); }

As illustrated below, the accuracy of the result may not match that obtained using the spectrum technique (though this could be subjective). Nonetheless, the approximation is still quite good and the function is approximately 66 times faster to compute than the spectrum version. For real-time applications, this equation, combined with the Stefan-Boltzmann law for assessing the total intensity emitted by the blackbody, might be more advantageous.